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ABSTRACT 

The dayside of HD 149026b is near the edge of detectability by the Spitzer Space Telescope. We report on 
eleven secondary-eclipse events at 3.6, 4.5, 3 x 5.8, 4 x 8.0, and 2x16 |xm plus three primary-transit events at 
8.0 [a.m. The eclipse depths from jointly-fit models at each wavelength are 0.040 ± 0.003% at 3.6 i,im, 0.034 
± 0.006% at 4.5 yim, 0.044 ± 0.010% at 5.8 ^m, 0.052 ± 0.006% at 8.0 |xm, and 0.085 ± 0.032% at 16 |xm. 
Multiple observations at the longer wavelengths improved eclipse-depth signal-to-noise ratios by up to a factor 
of two and improved estimates of the planet-to-star radius ratio {Rp/R^, = 0.0518 ± 0.0006). We also identify 
no significant deviations from a circular orbit and, using this model, report an improved period of 2.8758916 
± 0.0000014 days. Chemical-equilibrium models find no indication of a temperature inversion in the dayside 
atmosphere of HD 149026b. Our best-fit model favors large amounts of CO and CO2, moderate heat redistribu- 
tion (/ = 0.5), and a strongly enhanced metallicity. These analyses use BiLinearly-Interpolated Subpixel Sen- 
sitivity (BLISS) mapping, a new technique to model two position-dependent systematics (intrapixel variability 
and pixelation) by mapping the pixel surface at high resolution. BLISS mapping outperforms previous methods 
in both speed and goodness of fit. We also present an orthogonalization technique for linearly-correlated pa- 
rameters that accelerates the convergence of Markov chains that employ the Metropolis random walk sampler. 
The electronic supplement contains light-curve files and supplementary figures. 
Subject headings: planetary systems — stars: individual: HD 149026 — techniques: photometric 



1. INTRODUCTION 

Discovered in 2005 using Doppler measurements, the 
Saturn-sized extrasolar planet HD 149026b orbits (in 2.876 
days) a GOIV star that is larger (1.45 solar radii), and hot- 
ter (6150 ± 50 K) than most stars known to host transiting 
exoplanets. The planet's small radius and high average den- 
sity suggest that between 50% an d 90% of th e planet's mass 
must be in its rocky or icy core (iKnutson et al. 2009, here- 
after K09). Invoking current theories, it is difficult to form 
this exoplanetby gravitational instability (Sato et al. 2005). 

Shortly after detection, iFortney et al.l ([2006) computed 
models of the atmospheric temperature structure and spectra 
of HD 149026b. They suggested that the planet was a strong 
candidate for having a day-side atmospheric temperature in- 
version. The highly irradiated planet is hot enough to have 
gaseous TiO and VO molecules in the dayside atmosphere. 
These molecules are strong optical absorbers and had been 
previously s hown to cause temp erature inversions in model 
atmospheres (iHubeny et al.ll2003h . 

Beginning in 2005, we used the photometric channels of the 
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Spitzer Space Telescope (I Werner et al.l 12004') to observe HD 
149026b during secondary eclipse, when the planet passes 
behind its parent star, to cha r acteriz e the planet's dayside at- 
mosphere. iHarrington et al.l (120071 hereafter H07) found an 
8.0 \xm echpse depth of 0.084%+o o,'9, indicating the hottest 
brightness temperature (2300 ± 200 K) observed at that time. 
This temperature matches an insta ntaneous re-emission model 
(zero albedo) from Fortn ev et al.l (|2006) that exhibits a tem- 
perature inversion (which tends to enhance the planet/star 
contrast at 8.0 ^.m), thus suggesting the presence of absorbers, 
such as TiO and VO gas molecules, in the atmosphere. 

Charbonneau et al. (2006) observed two primary transits, 
when the planet passes in front of its parent star, using the 
Fred Lawrence Whippl e Observato r y tele scope through the 
Sloan g and r filters. IWinn et aP (120081) reported on five 
ground-based transits through Stromgren b and y fil ters at 
the Fairborn Observatory. In August of 2007, Nutzmai i et al.l 
(2009) used Spitzer to monitor a transit of HD 149026b at 8.0 
l^m. Carter et al. (2009) used the NICMOS detector on board 
the Hubble Space Telescope to observe four transits of HD 
149026b at 1.4 p.m. Their data have the best photometric pre- 
cision to date and, after combining their data with all previous 
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transit measurements, provide improved estimates of orbital 
parameters, mass, and radius. 

In 2008, K09 monitored the system for just over half an 
orbit to characterize the planet's phase variation at 8.0 |J.m. 
Their observations began slightly before the primary transit 
and finished slightly after the secondary eclipse. Using the 
final 7.2 hours of data, K09 report an eclipse depth of 0.041 1 
± 0.0076%, half that of H07. As part of their paper, K09 
reanalyzed the 2005 secondary-eclipse data and found eclipse 
depths ranging from 0.05 - 0.09%, though the lower values are 
preferred in most of their models. This large range of eclipse 
depths depends on the choice of systematic error model, fit- 
ting routines, and bad-pixel trimming methods. 

HD 149026b is an interesting planet given its extremely un- 
usual bulk abundances; the majority of the planet's mass must 
be in heavy elements, making the planet perhaps more akin to 
Uranus and Neptune than Jupiter and Saturn. In the solar sys- 
tem, a bulk composition that is enhanced in metals goes hand 
in hand with an atm ospheric composition enhanced in metals 
jMarlev etalj|2007h . This suggests that HD 149026b could 
have an atmospheric metallicity far greater than that of most 
transiting exoplanets. Verifying this by measurement would 
let us understand the makeup of this planet and the role of at- 
mospheric composition in determining temperature structure. 

In this paper we present Spitzer Space Telescope secondary- 
eclipse observations of HD 149026b that resolve the disagree- 
ment in eclipse depths at 8.0 [im, characterize the planet's 
dayside atmosphere, and further constrain its orbital and phys- 
ical parameters. We give detailed descriptions of our tech- 
niques and results because how one handles Spitzer 's system- 
atics can lead to best-fit parameters that disagree by more than 
1(7, as demonstrated in Section|5] 

Below, we describe the observations and data analysis, 
present a new method for modeling one of Spitzer' s system- 
atics, explain how we arrived at the final fits and compare the 
results to previously published work, discuss implications for 
the planetary emission spectrum and planetary composition, 
give improved constraints on the orbital parameters, and state 
our conclusions. 

2. OBSERVATIONS AND DATA ANALYSIS 

2.1. Observations 

We observed secondary eclipses of HD 149026b at 3.6, 4.5, 
5.8, a nd 8.0 |J.m with the Infrared Array Camera dpazio et alJ 
l2004l IRAQ and a t 16 |j.m using the Infrared Spectrograph's 
jHouck et al.ll2004l IRS) photometric blue peak-up array. The 
program also observed a primary transit at 8.0 |a.m. Including 
the four previously analyzed data sets labeled in Table [U we 
present fourteen observations spanning more than 3.5 years. 

2.2. POET Pipeline 

Our Photometry for Orbits, Eclipses and Transits (POET) 
pipeline produces systematics-corrected light curves using 
Sp/fzer-supplied Basic Calibrated Data, fits a multitude of 
models with a wide range of analytic forms for systematic ef- 
fects, chooses the best-fit model, and assesses the uncertainty 
of each free parameter. Below, we describe each of these steps 
in detail. 

We calculate the Julian date of each image at mid-exposure 
using the UTCS_OBS and FRAMT IME keywords in th e 
5;?/fzer-supplied headers. Following lEastman et al.l (120101) . 
we convert dates to Barycentric Julian Dates in the Coordi- 
nated Universal Time standard (BJDyjf^) using the JPL Hori- 



zons systerrQ to interpolate Spitzer' s position relative to our 
solar system's barycenter. Additionally, converting from UTC 
to the Barycentric Dynamical Time (TDB) standard addresses 
any discontinuities due to leap seconds. This paper reports 
both time standards to facilitate comparisons with previous 
work, which mostly does not apply the leap-second correc- 
tion, and to ease the transition to the more-accurate standard. 

The POET pipeline flags bad pixels (from energetic parti- 
cle hits and other causes) by grouping sets of 64 frames and 
doing a two-iteration, Aa rejection at each pixel location in 
the set. Stellar centers for photometry come from a Gaus- 
sian fit and 5 x interpolated aperture photometry (H07 Sup- 
plementary Information, SI) produces the light curves. We 
test a broad range of aperture sizes in 0.25-pixel increments 
and omit frames with bad pixels in the photometry aperture. 
The background, subtracted before photometry, is an average 
of the good pixels within the specified annulus centered on 
the star in each frame. 

2.3. S\i\tz&TL Systematics 

Exoplanet characterization require s photometric sta bility 
well beyond Spitzer's design criteria (F azio et 311120041) . De- 
tector sensitivity models vary by channel and can have both 
temporal (detector ramp) and spatial (intrapixel variability) 
components. The main systemati c effect at 3.6 and 4 .5 M.m 
is intrapixel sensitivity variations (ICharbonneau et an r2005). 
in which the photometry depends on the precise location of 
the stellar center within its pixel. We fit this systematic using 
the new BiLinearly-Interpolated Subpixel Sensitivity (BLISS) 
mapping technique described in Section [3] This technique 
maps the spatial sensitivity variations at high resolution. 

The 5.8, 8.0, and 16 |J.m arrays primarily suffer from tem- 
poral variability, attributed to charge trapping (K09) in the 8.0 
|j.m case. Weak spatial dependencies can also occur at these 
wavelengths (Stevenson et al. 2010), so we consider both sys- 
tematics when determining our best-fit model. Typically, we 
omit the initial portion of each light curve from the model fit 
to avoid the worst of the ramp effect and to allow the telescope 
pointing and detector to stabilize. The clipping parameter, q, 
defines the number of unmodeled points from the start of a 
data set. 

2.4. Pixelation 

Pixelation is an infrequently discussed systematic error in- 
herent to all array detectors. Sufficiently small stellar center 
motions between frames will not add or subtract pixels from 
an aperture, but these motions will cause the total flux within 
the aperture to vary. This means there is a (potentially large) 
range of stellar centers that utilizes the same set of aperture 
pixels, introducing a position dependence to the photometric 
sensitivity. We provide an illustrative example in Figure [T] 
and display the magnitude of this effect in Figure |2l A flux- 
conserving, subpixel image interpolation, combined with pre- 
cise centering and applied before photometry, mitigates the 
pixelation effect by decreasing its range and amplitude. As 
demonstrated by our BLISS maps in Section[3] uninterpolated 
photometry exhibits strong sensitivity peaks at one-pixel in- 
crements, while 2 X - and 5 x -interpolated pixels exhibit pro- 
gressively weaker peaks at 0.5 and 0.2 pixel increments, re- 
spectively, in each spatial direction. 

Pixelation is most apparent with small apertures placed 
on under-resolved point-response functions (PRFs) such as 

' http://ssd.jpl. nasa.gov/?horizons 
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TABLE 1 
Observation Information 



Label'' 


Observation Date 


Duration 


Frame Time 


Total Frames 


Spitzer 


Wavelength 


Previous 






[minutes] 


[seconds] 




Pipeline 


[|j.m] 


Publications'' 


HD149bs41 


August 24, 2005 


330 


0.4 


44352 


SIS 


.7.0 


8.0 


H07 


HD149bs51 


August 4, 2007 


386 


14 


1050 


SIS 


.18.0 


16 




HD149bs31 


August 13, 2007 


386 


0.4 


54080 


S18.7.0 


5.8 




HD149bp41 


August 14, 2007 


478 


0.4 


67008 


SIS 


.7.0 


8.0 


N09 & C09 


HD149bs52 


August 30, 2007 


386 


14 


1050 


SIS 


.18.0 


16 




HD149bp42 


Sept. 12, 2007 


386 


0.4 


54080 


SIS 


.7.0 


8.0 




HD149bsll 


March 10, 2008 


386 


0.4 


54080 


SIS 


.7.0 


3.6 




HD149bs42 


April 11,2008 


386 


0.4 


54080 


SIS 


.7.0 


8.0 




HD149bs21 


May 9, 2008 


386 


0.4 


54080 


SIS 


.7.0 


4.5 




HD149bp43 


May 11,2008 


499 


0.4 


70000 


SIS 


.7.0 


8.0 


K09 


HD149bs43 


May 12, 2008 


432 


0.4 


60500 


SIS 


.7.0 


8.0 


K09 


HD149bs32 


June 16, 2008 


386 


0.4 


54080 


SIS 


.18.0 


5.8 




HD149bs33 


March 13, 2009 


386 


0.4 


54080 


SIS 


.18.0 


5.8 




HD149bs44 


March 22, 2009 


386 


0.4 


54080 


S18.18.0 


8.0 





''HD149b designates the planet, p/s specifies primary transit or secondary eclipse, and ## identifies the wavelength and 
observation number. 

''H07 = IHarrington"etaI] <2007h , N09 = INutzman elal] {20091) , C09 = ICarteretaI] {20091) , and K09 = IKnutson et"al] 
(.20091) . 
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Fig . 1 . — Illustrative example of the pixelation elfect, which arises when 
small changes in stellar centers do not add or subtract pixels within the aper- 
ture but do change the collected flux. The left panel uses a solid blue circle to 
depict a photometry aperture centered at (15.5, 15.5), which is on the corner 
of four pixels (defined by dashed black lines). The 12 shaded pixels have 
centers (green dots) that fall within the aperture and are summed to deter- 
mine the stellar flux in this image. All of the incoming photons that land 
within the dashed blue circle (where there is a greater density of incoming 
photons) count towards the flux. So long as the stellar center falls within the 
small black circle, the photometry aperture will encompass the same green 
dots; thus, the specific pixels that contribute to the total flux will not change. 
One such example is depicted in the right panel, where we apply an offset 
of (0.15, 0.15) pixels. In this case, not all of the photons that fall within the 
dashed red circle count towards the flux. Instead, the shaded pixels count 
additional flux from photons that land outside the aperture (solid red circle) 
where the density of photons is less. The net effect is that, due to a change 
in centering and a non-uniform photon density, the shaded pixels in the right 
panel will record fewer incoming photons. The result is a position-dependent 
systematic called pixelation. 



Spitzer's. It may not be apparent in other situations, such as 
when the aperture contains almost all of the integrated PRF, 
when the centroid wander is small relative to the subpixel size, 
and when other noise sources dominate (e.g., systematics and 
variable PRFs). Increasing the aperture size lessens the pix- 
elation effect by decreasing the fraction of uncaptured light 
outside the aperture, but doing so may decrease the signal- 
to-noise ratio (S/N) by increasing the amount of background 
noise included in the aperture. Thus, choosing the best aper- 
ture size may introduce this position-dependent systematic. 
We have determined the intrapixel variability at 5.8, 8.0, and 
16 i^m to be a pixelation effect, previously reported at 5.8 )i.m 
by .Stevenson et al. (2010.) and at 8.0 \im by .Anderson et al.l 



• • Binned Flux 
— BUSS Map 




Pixel Postion in y 



Pixel Postion in x 



Fig. 2. — Projected flux from HD149bs31 integrated along the x (left) and y 
(right) axes. The non-uniform flux in both panels is clear evidence of pixela- 
tion. We use 5 X -interpolated aperture photometry, which results in 0.2-pixel 
spacing between peaks. In the left panel, low-order polynomial models would 
fit the pixelation effect poorly at the peaks; the BLISS map (see Section [3) 
has no such limitation. The systematic is weakly constrained near the edges 
due to low sampling, as indicated by the large error bars. Whether a specific 
point on a pixel is a local maximum or minimum (due to pixelation) is a func- 
tion of aperture size, which defines which subpixels to include in the aperture 
at any given point on the detector. 

(|20T1 . 

There are several ways to correct pixelation. First, one 
could shift model PRFs to match each frame's precisely de- 
termined stellar center Dividing the stellar flux by the PRF 
flux in the aperture should remove the effect, but requires a 
highly accurate model PRF. Second, using smaller subpixels 
could decrease the amplitude of the pixelation effect until it is 
insignificant relative to the noise, but there is a limit: interpo- 
lation can only approximate the information destroyed when 
photons fall into the detector's finite-sized pixels. Third, if 
the pointing is sufficiently consistent and compact, one could 
choose an interpolation factor that happens to place the flat 
portion of the pixelation response on the stellar centers (since 
the peaks in Figure |2] move with different interpolation fac- 
tors). Fourth, with high-precision centering, one could use a 
series of images taken at slightly different positions to model 
the position sensitivity analytically or with pixel-mapping 
techniques such as BLISS, but one would first have to re- 
move any time-dependent components from the light curve 
model (see Section 12.51) . The accuracy of these models de- 
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pends on the centering and photometric precisions, the former 
being ^0.01 pixels for 0.4 secon d IRAC subarray expo sures 
of bright sources (see below and'Steve nson et aTIl2010l) . We 
test for pixelation in all data sets using BLISS mapping, which 
corrects the effect when it is significant. 

2.5. Light-curve Modeling 
The full Ught-curve model is: 

F(x,y,t) = FMt)R{t)M{x,y)V(v)P(p), (1) 

where F(x,y,t) is the measured flux centered at position {x^y) 
on the detector at time t, fs is the (constant) system flux 
outside of secondary eclipse or primary transit, E{t) is the 
primary-transit or secondary-eclipse model component, R{t) 
is the time-dependent ramp model component, Mix^y) is the 
position-dependent intrapixel model component or sensitivity 
map, V(v) is the visit sensitivity as a function of visit frame 
number v, and P{p) is the flat-field correction at position p. 
Below, we discuss some of these components in more detail. 

2.5. 1 . Eclipse and Transit Models 

The uniform - source and small-planet equations from 
iMandel & Agoll (I2002I) describe the secondary-eclipse and 
primary-transit model components, E{t), respectively. Tran- 
sit light curves at 8.0 |J.m exhibit weak limb darkening that is 
not well constrained by fitting li mb-darkening mode ls to the 
data. We foUow the method of iBeaulieu et alJ (l2008i) in de- 
riving limb-darkening coefficients for HD 149026. Spitzer's 
8.0 |xm spectral response curve weights t he intensities of a 
Kuru cz ATLAS stellar atmosphere model jCastelli & Kurucj 
I2OOI Teff = 6250 K, logig) = 4.5 cgs, and [M/H] = 0.3), given 
as a function of wavelength and angle from the star's center A 
least-squares minimization of the resulting curve determines 
the non-linear limb-darkening coefficients (Claret 2000, 01-04 
= 0.51477, -0.80525, 0.75683, -2.6168), which are then fixed 
for the three transit light-curve fits. 

2.5.2. Ramp Models 

We consider a multitude of ramp equations, R{t), all of 
which stem from three basic forms: exponential, logarithmic. 



and/or polynomial. 

R(t) = l±e-''°'^''' (2) 

R{t)=l±e-">'^'''+r2(t-0.5) (3) 

R(t)= l±e-">'^"+r2{t -0.5) + r^it-O.Sf (4) 

R(t)=l±e-'''''^'''±e-''"^''' (5) 

/;(f)= l + n(f-0.5) + r6ln(f-fo) (6) 

7^(0= l + ri(f-0.5) + r2(f-0.5)^ + r6ln(f-fo) (7) 

R(t)=l + r6lnit-to) + n[lnit-to)f (8) 

R(t) = 1 + ri(f - 0.5) + r2(t - 0.5 f + ln(f - to) + n[ln(t - to)f 

(9) 

R(t) =1 + nit -0.5) (10) 
R(t) =1 + nit -0.5) + r2it -0.5)^ (11) 



The time, f, is in units of phase (days) for secondary-eclipse 
(primary-transit) events. We use "H-" and "-" subscripts in Eqs. 
|2]-|5]to denote the corresponding functional form. For exam- 
ple, Eq. ^ describes an exponentially decreasing, asymptot- 
ically constant ramp while Eq. |2]_ describes an exponentially 
increasing, asymptotically constant ramp. 

There is a physic al interpretation applicable to the rising 
exponential ramps ( Agol et al]|2010 ). Consider a population 
of charge traps due to an impurity in the detector's infrared 
material and a flux of photoelectrons through the material. 
The traps collect some fraction of electrons, releasing them 
randomly with some characteristic time scale that depends on 
the impurity. Bright sources saturate the traps, decreasing the 
fraction of captured electrons and raising the detected signal 
of a steady source according to the asymptotic rising exponen- 
tial function in Eq. |2]_. A double rising exponential (Eq. |5l_) 
approximates a rapidly saturating PSF core and slowly satu- 
rating wings (Knutson et al. 2008). It could also represent two 
impurities; the very short ramps of the HD149bs41 data set's 
visit sensitivity suggests the presence of an impurity that has 
a very short time scale and releases many of its electrons over 
the course of one cycle (^30 minutes, see H07 Supplemen- 
tary Figure 6). We tested for a common set of characteristic 
time scales in all of the rising exponential free parameters; 
however, even for the same planet in the same array, we did 
not achieve reasonable fits with all data sets. This may be due 
to inadvertent and inconsistent pre-flashing by the objects ob- 
served prior to our own observations. Despite its potential for 
providing a physical explanation for the ramps, the rising ex- 
ponential does not alwa ys pr ovide the best model according 
to the criteria in Section 12.61 This may be due to our fitting 
the final photometry and not the individual pixels' data and/or 
to the pointing instability producing unsteady illumination at 
the precision levels relevant here. 

Agol et al. (2010) advocate using a double rising exponen- 
tial (Eq.|5]_) for 8.0 |xm data due to its improved fit and smaller 
residuals, weaker dependence on aperture size , and less sen- 
sitivity in the clipping parameter Similarly, iKnutson et al.l 
( I2OI ll) find that Eq. |2]_ is sufficient for their 8.0 |i.m pre- 
flashed data sets, while Eq. |5]_ is necessary for their non- 
preflashed data. We test all relevant ramp equations on data 
sets that exhibit time-dependent systematics and orthogonal- 
ize any correla ted para meters that inhibit convergence (see 
Sections 12.5.31 and 12.7b . Upon doing so, we find that we 
cannot corroborate the claims by lAgol et all (120101) . Equa- 
tion |5]_ should not be used for all 8.0 |^m data sets because 
we typically find correlations between the eclipse depth and 
its ramps parameters that can double the latter' s uncertainty 
relative to other models (see HD149bp42). Rather, we rec- 
ommend a comprehensive examination of all relevant ramp 
equations before selecting the final model. See Sections |4] 
and|5]for discussion relevant to particular events. 

2.5.3. Orthogonalization 

The exponential model components have one difficulty: 
the Markov-Chain method used to assess the uncertainties 
does not converge to the posterior probability distribution, 
according to the criteria discussed in Section 12.71 even af- 
ter tens of millions of iterations. The problem is a correla- 
tion between the expone ntial ramp parameters. Several re- 
cent a nalyses (e.g., K09. IStevenson et £11120 lOt TCampo et al.l 
I2OI ll) have "solved" the problem by fixing one of the ex- 
ponential parameters. This is effectively profiling (slicing) 
rather than marginalizing (integrating) over the posterior pa- 
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rameter distribution, an approach that ignores co rrelated er- 
rors a nd may reduce the calculated error bar (see iPress et aTl 
120071 Figure 15.6.3 and related text). In one case discussed 
below (HD149bp42), fixing parameters incorrectly decreased 
the calculated eclipse depth uncertainty by 50%. There are 
two legitimate escapes from the problem. The first is to write 
a Monte Carlo sampler that explores the phase space more 
intelligently than the Metropolis random walk we and many 
others use. The second is to re-cast the equations in a form 
that eliminates the non-linear correlation of the ramp parame- 
ters. 

For this paper, we have expediently chosen the second ap- 
proach. The version of Eq.|2]_ presented here produces a more 
linear correlation between the ro and ri parameters than the 
original version in H07, whose exponent was -r(t - to). In 
some cases this converges quickly and the job is done. In oth- 
ers it still does not converge, so we run at least 10^ iterations 
to sample the posterior distribution and then rewrite the model 
with a change of variables that orthogonalizes the most cor- 
related parameters. This method does not modify the number 
of free parameters, nor involve any interpolations or approx- 
imations. For the selected correlated parameters, the orthog- 
onalization shifts the origin to the center of the distribution, 
divides each parameter by its standard deviation to give a uni- 
form scale in all directions, and rotates the subspace to min- 
imize correlations (see Figure |3]l. A routine that prepares for 
principal components analysis (PCA) finds the transformation 
matrix from the distribution sample (we do not actually per- 
form the PCA). We rerun the Markov chain using the new 
equations until it converges according to the criteria given be- 
low. Then we transform the points back to the familiar equa- 
tions to assess parameter uncertainties. A simple example of 
a two-parameter orthogonalization of Eq. |2]is: 

R(^t)= 1 zbe"*^"'' cos(e)-sm(e)] ^ci [( sin(e)+cos(6l)] ^ 2) 

where cq and ci are the ramp parameters in the rotated frame 
and 9 is the rotation angle between coordinate systems. The 
arctangent of the slope of the best-fit line through the initial 
sample, projected into the r^-ri plane, determines 9. Such 
approximate orthogonalization of the posterior distribution 
has long been st andard practice for improvin g MCMC per- 
formance (e.g.. Hills & Smith 1992; Brooks 1998). Similar 
dis cussion can be found in Connollv et al. (1995), Pal (2008). 
and lCowan et"!!] ( I2009I) . 

In effect, this method is the same as the first, accom- 
plished by rotating the data and using the original sampler 
rather than writing a new one. This method works best 
for linearly-correlated parameters and achieves moderate suc- 
cess with more exotic correlations. Converting to curvi- 
linear coordinates or implementing manifold learning algo- 
rithms from nonlinear dimensionality reduction may offer fur- 
ther improvements in convergence time for the extreme cases 
£ee& Verleysen 2007|). 

2.5.4. Flat Field (Position) Sensitivity Models 

Most of the data sets presented here follow the standard 
time-series observing practice of keeping the object fixed 
to one location on the array (staring) to minimize position- 
dependent sensitivity effects. However, the H07 observa- 
tion (HD149bs41) cycled through nine different nod positions 
(p = - 8) in an attempt to use the unobserved positions in 
each frame to make a high-quality flat field that would cor- 
rect the entire data set. This approach was unsuccessful and 
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Fig. 3. — Two-parameter orthogonalization example for HD149bs21 with 
histograms. The physical parameters (top panels) show a strong, non-linear 
correlation and asymmetric histograms; however, the orthogonalized parame- 
ters (bottom panels) are nearly uncoirelated and have symmetric histograms. 
Running a Markov-chain method with the orthogonalized parameters reduces 
the convergence time. 

not repeated. Each position in the HD149bs41 data requires 
a flat-field correction, P(p), to account for the difference in 
pixel sensitivity. To eliminate correlations with the system 
flux (i.e., keep the correction from floating), we require the 
mean of all of the corrections to equal unity. We do this by 
freely varying p = 1 - 8 and equating /:> = to the number of 
positions minus the sum of the other corrections. 

For both 16 i.im data sets, the telescope reacquires the target 
at one third and two thirds of the way into the observing runs. 
This action is similar to a nod because after reacquisition, 
the three sets of measured stellar centers are non-overlapping. 
We apply the same model component to these data as with 
HD149bs41, but only use three flat-field correction parame- 
ters (/? = 0,1,2). 

2.5.5. Visit Sensitivity Model 

With HD149bs41, Spitzer completed twelve cycles through 
the nine nod positions mentioned above for a total of 108 vis- 
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its. Each visit has a briefer and steeper ramp compared to the 
overall ramp, R{t). As discussed in their SI, H07 use a 12- 
knot spline to model the visit-sensitivity effect, V{v). They 
fix the final three knots to unity and allow the remaining nine 
parameters to vary. We model the visit sensitivity ramp using: 

y(u) = ui •ln(u-uo)+l, (13) 

which is identical in form to the model component used by 
K09. The only difference is that K09 use time as the indepen- 
dent variable while we use visit frame number, v. In either 
case, the independent variable resets to zero upon moving to 
a new nod position. 

2.6. Model Selection 

For each data set, we test different photometry aperture 
sizes, detector ramp model components, and intrapixel sen- 
sitivity model components looking for the best combination. 
One must be careful in assessing which model is the "best" 
because x^-like comparisons must derive from the same data 
set and different photometric apertures produce different data 
sets for this purpose. We use the standard deviation of the 
normalized residuals (SDNR) when comparing models of the 
same analytic form to different data sets. Once we have iden- 
tified the best ap erture size, we use the Bayesian Information 
Criterion (BIC, Schwarz" 1978t lLiddlell2007t IStevenson etaTI 
12010: Campo et al. 2011J to compare models with different 
numbers of free parameters: 

,2 

BIC= -^ + klnn. (14) 

Here, e, and ct,- are the residual and uncertainty of the /* data 
point, k is the total number of free parameters, and n is the 
number of data points used in the model fit. The Spitzer- 
supplied uncertainty frames are typically overestimated, so 
we scale cr, such that the reduced xi, equals unity for 
our best-fit model. Using the unsealed ct,- values in a least- 
squares minimizer improperly weights the data and results in 
a sub-optimal fit. When selecting between competing models 
for an observed dataset, what matters is how the models com- 
pare in predicting the actually observed data, not how variable 
individual model fits may be among hypothetical realizations 
of the data. The BIC is defined so that the marginal likeli- 
hood ratio — the ratio of predictive probabilities for the ob- 
served data — is approximately gO.SAB/c jsjote that what mat- 
ters is the dijference in BIC values (and thus the difference of 
maximum likelihood or minimum values), not the absolute 
values. Information criteria such as BIC may not apply when 
comparing fits with intrapixel maps to those with polynomial 
model components. See Appendix lAl for more discussion on 
the matter. 

2.7. Error Estimation 

We explore phase space and estimate uncertainties using a 
Mai-kov-chain Monte Carlo (MCMC) routi ne follow i ng the 
Metropolis random walk algorithm. See Ca mpo et al.l (1201 lb 
for more discussion on MCMC. The POET pipeline can test 
any combination of systematic model components, comput- 
ing the SDNR and BIC for each fit, and can model multi- 
ple events simultaneously while sharing parameters such as 
the eclipse midpoint and depth. Each MCMC run begins 
with a least-squares minimization, a rescaling of the Spitzer- 
suppHed uncertainties, a second least-squares minimization 



using the new uncertainties and, finally, at least 10^' MCMC 
iterations to populate the posterior parameter distributions, 
from which we derive parameter uncertainties. We study pa- 
rameter correlation plots and the posterior distribution to en- 
sure that we have a reliable result, then publish them so that 
others may evaluate our work and compare to their own. We 
tes t for convergence every 10^ steps, terminating only when 
the lGelman & Rubinl(ll992 !) diagnostic for all free parameters 
has dropped to within 1 % of unity using all four quarters of 
the chain. We also examine trace and autocorrelation plots of 
each parameter to confirm conv ergence visually . We estimate 
the effective sample size (ESS. IKass et al.lll998h and autocor- 
relation time, T, for the /* free parameter as follows: 



*:=I 

where m is the length of the MCMC chain, pk{Oi) is the au- 
tocorrelation of lag k for the free parameter 6',, and k' is the 
cutoff point such that < 0.01. We use the longest auto- 
correlation time from each event to determine the number of 
steps between effectively independent samples for thinning 
each MCMC chain. The distribution of thinned samples quan- 
tifies parameter uncertainties. 

3. BLISS MAPPING TECHNIQUE 
3.1. Background 

The change in pixel sensitivity with respect to stel- 
lar position on the detector is a well known systematic 
with the Spitz er Space Telescope dCharbonneau et all 120051 : 
iKnutson et al.l [2008). This effect is particularly strong in 
IRAC's 3.6 and 4.5 M.m channels but has also been seen at 5.8, 
8.0, and 16 |J.m (Stevenson et al. 2010; Anderson et al. 201 1). 
The position sensitivity in the latter channels is due to a pix- 
elation effect (see Section |2] for a description) rather than an 
actual change in sensitivity over the pixel surface. In this sec- 
tion, we present a new technique for modeling these position- 
dependent systematics, called BiLinearly-Interpolated Sub- 
pixel Sensitivity (BLISS) mapping. 

The most common method for removing the intrapixel 
variability is to fit a polynomial in both spatial directions 
(IKnutson et al.ll2008l) . The polynomial order typically ranges 
from quadratic to sextic and may include cross terms. Other 
variations of this modeling method have applied multiple 
polynomials, one for each pixel quadrant, in order to find the 
best fit. Polynomial methods work reasonably well for data 
sets with small stellar position wander, resulting in a smoothly 
varying intrapixel sensitivity. An analysis becomes exceed- 
ingly complicated if the variation is not smooth or if strong 
correlations arise between parameters. These complications 
can increase uncertainty estimates in the best case scenario 
and, in the worst case, lead to inco rrect results. 

A new approach, pioneered by iBallard et alj (|2010, here- 
after BIO), attempts to map the intrapixel variability on a 
subpixel-scale grid without assuming a specific functional 
form. For their particular light curve, they bin their flux and 
stellar positions into 20-second bins (^ 145 points/bin) before 
computing a sensitivity correction for each binned point. Each 
correction considers a set of flux values that does not include 
in-transit frames or frames from the current binned position. 
This set of flux values is Gaussian- weighted in both spatial di- 
rections relative to the position of the binned point being cor- 
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reeled, summed to a single value, then normalized by dividing 
by the summed Gaussian weighting function. This method ef- 
fectively models not only the large-scale features in their data, 
but also a smaller-scale "corrugation" effect that a low-order 
polynomial cannot remove. 

Moving away from a polynomial model is an excellent con- 
cept; however, this particular implementation has some draw- 
backs that limit its scope. For instance, ignoring the points 
during secondary eclipse requires that the out-of-eclipse por- 
tion of the data set be significantly longer than the in-eclipse 
portion, which is atypical in primary-transit and secondary- 
eclipse observations. Also, the weighting function computes 
too slowly to be used in a MCMC routine and is even slow 
when using a minimizer The calculation would be even 
slower if the data were not binned into relatively long, 20- 
second time intervals. Figure |4] shows 120 seconds of verti- 
cal pixel positions from HD149bsll, which has a 0.4-second 
exposure duration vs. 0.1 seconds for BIO. The stellar cen- 
ter can vary significantly over a 20-second interval, indicat- 
ing that positional information is lost when binning over such 
time scales. 
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> 
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Fig. 4. — Vertical pixel positions from HD149bsll. We can track the mo- 
tion of the stellar center to high precision (~0.01 pixels) over a duration of 
only a few seconds. These oscillations are much faster than Spitzer's 1/2-1 
hour oscillations, initially reported by Charbonneau et al. (2005). In a span 
of 20 seconds, the stellar centers can vary by > 0.1 pixels. All positions are 
zero-based. 



The idea of using a spline to interpolate position-dependent 
systematics stems from observed fine-scale sensitivity varia- 
tions in some of our data sets that cannot be modeled by a low- 
order polynomial (see the pixelation effect in Figures |2] and |5] 
for examples). We initially attempted to map the pixel surface 
using a bicubic spline because we wanted a smoothly-varying 
model; however, this type of interpolation is prohibitively 
computationally expensive. A typical secondary-eclipse ob- 
servation spans a 0.3 x 0.3 pixel region. Placing knots at 
0.05-pixel intervals requires 49 free parameters. Polynomial 
models typically require less than 10 parameters. Adequately 
describing the fine-scale sensitivity variations requires a large 
number of knots, but varying all of these knot parameters at 
each step of an MCMC routine leads to extremely slow con- 
vergence. Our new BLISS mapping technique circumvents 
the problem of slow convergence by directly computing the 
knot parameters, rather then allowing them to vary freely. 
Thus, we can use >1000 knots to map the pixel surface at 
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Fig. 5. — BLISS maps illustrating the position-dependent pixelation effect. 
The maps use 1 X - (left), 2x - (center), and 5 X - (light) interpolated, 2.5 pixel- 
aperture photometry. Redder (bluer) colors indicate more (less) flux within 
the aperture. The bin size is 0.01 pixels for all maps. The horizontal and 
vertical black lines depict pixel boundaries. Without subpixel interpolation, 
the pixelation effect is significant, but it is progressively reduced with 2 X - 
and 5 X -interpolated photometry. For time-series data such as these, one can 
calculate a BLISS map to correct for pixelation. 

high resolution (see Figure|5]l. 

3.2. Implementation 

In seeking methods that compute faster than bicubic in- 
terpolation, we give up the constraint of differentiability. 
Nearest-neighbor interpolation (NNI) is simple, assigning 
each point the value of its nearest knot. Bilinear interpola- 
tion (BLI) is a straightforward calculation (see Eq. [TSI l that 
maintains continuity over the pixel surface, unlike NNI, and, 
given sufficiently precise centering, should be more accurate. 
Using BLI, the flux at a point (x^y) is: 



M{x,y)=Fip{xx,yx) 
+ Fw(x2,y\) 



+ Fip(xi,y2) 



+ FiAx2:yi) 



)— 
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y\ 



(16) 



X2-XI y2-y\ 

This is a distance-weighted average of the flux of the four 
nearest knots, Fi^ixi^yj), where / and j are horizontal and ver- 
tical indices for a rectangular grid of knots. This method 
computes faster than bicubic interpolation and may achieve 
comparable smoothness within the errors with less computing 
time simply by increasing the number of knots (see Figure|2]). 

We create a rectangular grid of knots that spans the range of 
centers in x and y. Each point in the data set associates with 
its nearest knot. For BLI, we compute the distances from each 
point to its four nearest knots, for Eq.[T6] If one or more knots 
in Eq.[T6ldoes not have any assigned points, we use NNI there 
instead, or the calculation would fail. This usually only occurs 
near the boundary of the grid of knots. We precompute the 
knot associations and distances prior to initiating the MCMC 
as they remain constant from iteration to iteration. 

We do not treat the knots as MCMC jump parameters. 
Rather, we step all other free parameters from Eq. [T] gen- 
erate a new model using these new jump parameters, then 
divide the observed flux by the new model {F]p(x,y) = 



8 



Stevenson et al. 



F„hs/F,E{t)R{t)V{v)P{p)). Hypothetically, the residuals of 
Fwix^y) contain only position-dependent flux variations. The 
flux value of a particular knot is the mean of Fi^ix^y) for the 
points associated with that knot. We also tried median and 
weighted average knot values but the results did not improve 
and these calculations are much slower Next, we generate 
the sensitivity map (M[x,3'], Figure |6]l by interpolating the 
flux from the knots to all of the observed points using BLI 
and/or NNI. We tested various weighted smoothing functions 
when generating the sensitivity maps but, again, there was 
no improvement in the results. Finally, M{x,y) enters Eq. [1] 
for comparison with the observed flux to obtain an estimate 
of the goodness of fit and determine the MCMC acceptance 
probability. This process repeats for each step of the MCMC 
routine or minimizer. 
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Fig. 6. — BLISS map and pointing histogram of HD149bsll. Top: Red- 
der (bluer) colors indicate higher (lower) subpixel sensitivity. The horizontal 
and vertical black lines depict pixel boundaries. Bottom: Colors indicate the 
number of points in a given bin, which, in this case, is 0.015 pixels in length 
and width. By recalculating the map at each step of the MCMC or mini- 
mizer, this technique substantially improves on that of iBallard et ah >2010i) . 
and beats all tested functional fits. 



3.3. Determining The Optimal Bin Size 

The accuracy of the BLISS mapping technique depends 
critically on the bin size, or resolution in position space; how- 
ever, there is a trade-off between bin size and speed. Decreas- 
ing the bin size requires more knots and runs slower but may 
be necessary to adequately resolve sensitivity changes on the 
pixel surface. There is, however, a practical limit to how small 



the bins can be. A bin for every measurement will always pro- 
duce a perfect fit, resulting in a negative number of degrees 
of freedom and leaving the eclipse parameters unconstrained. 
Bin sizes must be small enough to resolve real, small-scale 
variations on the pixel surface but large enough to mix in- and 
out-of-eclipse points. This mixture helps to minimize corre- 
lations between the eclipse parameters and the knots in the 
sensitivity map. 

To establish the optimal bin size, we consider a range of 
bin sizes for both BLI and NNI using the 3.6 |J.m data set 
(HD149bsll), which has the strongest position-dependent 
systematic. We draw several conclusions from the results 
in Table |2l First, the SDNR (our measure of goodness of 
fit) decreases with decreasing bin size, indicating a better fit. 
Unfortunately, the SDNR decreases indefinitely, so it cannot 
constrain the minimum bin size. Second, BLI fits the data 
better than NNI for bin sizes greater than 0.015 pixels. The 
opposite is true for smaller bin sizes, which is counterintu- 
itive because BLI should always outperform NNI, assuming 
no uncertainty in the position. Thus, the bin size at which 
NNI outperforms BLI is indicative of the centering precision 
for a particular data set. We estimate the precision for our 
analysis of HD149bsl 1 to be 0.009 pixels in x and 0.007 pix- 
els in y by calculating the root-mean-squared (RMS) frame- 
to-frame position difference. This agrees well with Figure |4] 
and is consistent with the cross-over bin size of 0.015 pixels, 
where NNI outperforms BLI. Last, we place an upper limit on 
the bin size by noting that the eclipse depths become incon- 
sistent with each other for pixel sizes > 0.050 using BLI and 
> 0.020 using NNI. 

We conclude that, whenever possible, BLI should be used 
with a bin size that is independent of the eclipse depth and 
has a lower SDNR than NNI. We have found cases where 
BLI is never better than NNI. In those instances, the position 
dependence is so weak that the intrapixel model component 
is unnecessary. A better fit with BLI, compared to NNI, is 
thus a good indicator that a position-dependence systematic is 
present in a given data set. For the test cases shown in Table 
121 using BLI with a bin size between 0.015 and 0.020 pixels 
is recommended based on our criteria. 

Precise centering is important for this method because im- 
precision limits the small est meaningful bin siz e. Our prelim- 
inary work (see the SI of lStevenson et al.ll2OI0l) indicates that 
the Gaussian and least-asymmetry centering methods are bet- 
ter than the center of light; additional work is in preparation. 

3.4. Comparing Intrapixel Models 

To compare the BLISS mapping technique with other in- 
trapixel methods, we fit six different intrapixel models to the 
HD149bsll data set. These models are quadratic, cubic and 
sextic polynomials (including lower-order cross terms), B lO's 
new weighted sensitivity function (fixing CTj to 0.021 and <jy 
to 0.0079), BLI, and NNI. The eclipse depth in BlO's model 
is slightly shallower than the other models, which are all well 
within 1(7 of one another, but the uncertainties are essentially 
identical (see Table [3]). The BLISS models show significant 
improvement in SDNR compared to the others. We cannot 
use the BIC to compare the BLISS model components with 
the polynomial model components for the reasons discussed 
in AppendixlAl 

BLISS mapping represents a substantial improvements over 
polynomial model components because BLI and NNI can 
model real structure (such as pixelation) that cannot be mod- 
eled with low-order polynomials, they encounter fewer corre- 
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TABLE 2 

BLISS Map Test - Variable Bin Size 



Model 


Bin Size 


SDNR 


Eclipse Depth 




[Pixels] 




[%] 


BLI 


0.100 


0.0028736 


0.063 ± 0.003 


BLI 


0.050 


0.0028222 


0.043 ± 0.003 


BLI 


0.020 


0.0028076 


0.040 ± 0.003 


BLI 


0.015 


0.0028031 


0.040 ± 0.003 


BLI 


0.010 


0.0027917 


0.040 ± 0.003 


BLI 


0.005 


0.0027403 


0.040 ± 0.003 


BLI 


0.002 


0.0024796 


0.039 ± 0.003 


NNI 


0.100 


0.0029435 


0.071 ± 0.003 


NNI 


0.050 


0.0028527 


0.054 ± 0.003 


NNI 


0.020 


0.0028116 


0.044 ± 0.003 


NNI 


0.015 


0.0028024 


0.041 ± 0.003 


NNI 


0.010 


0.0027865 


0.041 ± 0.003 


NNI 


0.005 


0.0027109 


0.040 ± 0.003 


NNI 


0.002 


0.0023773 


0.039 ± 0.003 



TABLE 3 

BLISS Map Test - Comparing To Other Intrapixel Models 



Model 


Bin Size 


SDNR 


Eclipse Depth 








[%] 


Ballard 


2.4 seconds" 


0.0028230 


0.034 ± 0.003 


Cubic 




0.0028180 


0.039 ± 0.003 


Sextic 




0.0028157 


0.040 ± 0.003 


Quadratic 




0.0028186 


0.041 ± 0.003 


BLI 


0.015 pixels 


0.0028031 


0.040 ± 0.003 


NNI 


0.015 pixels 


0.0028024 


0.041 ± 0.003 



" Longer bin sizes were considered but produced worse results. Shorter 
bin sizes are prohibitively expensive to compute and are below the limit of 
detectable motion. 

lations between free parameters, and require fewer iterations 
to assess parameter uncertainties. 

4. PRIMARY TRANSIT FITS AND RESULTS 

We present the scaling of the RMS model residuals vs. bin 
size (a test of correlation in time) in Figure|7] the best-fit tran- 
sit light curves in Figure [8] a comparison between two fits in 
Table 21 and the full set of best-fit transit parameters in Ap- 
pendix|B] The electronic supplement contains light-curve files 
and supplementary figures. Below, we discuss each observa- 
tion to explain how we arrived at the final results. 

4. 1 . Three Fits at 8.0 fim 

At each 0.25-pixel increment in photometry aperture size, 
we model the time-dependent systematic with the ramps listed 
in Eqs.|2]-[II] An aperture size of 3.5 pixels produces the low- 
est SDNR values for all ramps and for all transit events. We 
estimate the background flux using an annulus from 7 to 15 
pixels centered on the star. We follow the method described 
in Section [33] when determining the optimal bin sizes of the 
BLISS maps. 

4.1.1. HD149bp41 

There are 26 consecutive frames (19494 to 19519) shifted 
horizontally by exactly one pixel; we flag these frames as bad. 
After clipping the first 5,000 data points (^33.3 minutes, q 
= 5,000), many of the ramps fit the time-dependent system- 
atic equally well, but the fits exhibit a large range of radius 
ratios (Rp/R^ = 0.0494 to 0.0517). Of the top three models 
shown in Table |4] Eq. |3]has the lowest BIC value and favors 
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Fig. 7. — RMS residual flux bin size for three HD 149026b transits. 
Blacli vertical lines at each bin size depict 1 ct uncertainties on the RMS resid- 
uals (RMS/ V^N, where N is the number of bins). The red line shows the 
predicted standard error for Gaussian noise, the dotted vertical blue line indi- 
cates the ingress/egress timescale, and the dashed vertical green line indicates 
the transit duration timescale. Any RMS residuals that are several a above 
the red line would indicate correlated noise at that bin size. When consider- 
ing the effects of congelations on transit depth, the bin size of interest is the 
transit duration and not the ingress/egress time. The shorthanded legend la- 
bels correspond to the last three characters in each event's label (e.g., p41 = 
HD149bp41). 

a moderate rad i us rat io of 0.0502 zfc O.OOl 1 . For comparison, 
iNutzman etall ( l2009h and lCarter etaTI (120091) report values of 
0.05158 ± 0.00077 and 0.5188 ± 0.00085, respectively To 
model the ramp, both adopt a quadratic function in ln(r) (Eq. 
[8]l with to fixed to a time a few minutes prior to the first obser- 
vation. Fixing parameters can cause a MCMC run to under- 
estimate uncertainties of the remaining free parameters (see 
Section |2. 5. 3b . Instead, we orthogonalize the correlated pa- 
rameters (system flux and all three ramp parameters in Eq. |3]l 
to provide a coordinate system in which our MCMC routine 
can sample efficiently. There is a weak position dependence 
(see Figure |8]l in both x and y directions that we model with 
the BLISS mapping technique using 0.030-pixel bins. 

4.1.2. HD149bp42 

For all ramps with reasonable fits, we find that the planet- 
to-star radius ratios range from 0.0444 to 0.513. The three 
best models appear in Table |4] with a range of uncertainties 
in Rp/Ri,. The ramp parameters from Eq. |5l_ correlate most 
strongly with the radius ratio, resulting in an uncertainty that 
is twice that of Eq. [8] Fixing ramp parameters in Eq. |5l_ er- 
roneously improves the radius ratio uncertainty to 0.008%. 
Equation [8] has the lowest BIC value, so we use it to fit the 
full data set iq = 0). The data exhibit only a minor posi- 
tion dependence in the x direction; however, the significant 
improvement in SDNR indicates that we should include the 
BLISS map during the joint model fit. 

4.1.3. HD149bp43 

Unlike the previous two data sets, HD149bp43 was pre- 
flashed to mitigate the ramp effect (see K09 for details). Thus, 
we do not need to clip a significant initial portion of the data 
set or use a double-rising exponential. Instead, we clip only 
the first 1000 data points (^6.7 minutes) and use Eq. |2]_ to 
model the time-dependent systematic. As seen in Table 21 
the HD149bp43 transit depth is consistently deeper than the 
other two data sets. The Rp/Ri, parameter is independent of 
our choice of q but is dependent on the choice of ramp model 
components, ranging from 0.0516 to 0.0536. BIC favors the 
deepest transit depth, resulting in a radius ratio that is larger 
than other best-fit ratios by ^ 4a. For comparison, K09 use 
Eq. in with a fixed to parameter and report a radius ratio of 
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Fig. 8. — Raw (left), binned (center) and systematics-corrected (right) primary-transit light curves of HD 149026b at 8.0 |xm. The results are normalized to 
the system flux and shifted vertically for ease of comparison. The colored lines are best-fit models, the black curves omit their transit model components, and 
the error bars are la uncerta intie s. The shorthanded legend labels correspond to the last three characters in each event's label (e.g., p41 = HD149bp41). The 
pixelation effect (see Section|T4) is most prevalent in HD149bp43. 



0.05253 ± 0.00076. We achieve the same underestimated un- 
certainty using a similarly constrained to parameter. A rela- 
tively strong position dependence is evident in Figure[8]and is 
modeled with a BLISS map using 0.050-pixel bins. 



TABLE 4 
Individual Transit Model Fits 



Label 


R(t) 


M(x,y) 


SDNR 


ABIC 




HD149bp41 




BLI 


0.0083449 


0.0 


0.0502 ± 0.0011 


HD149bp41 


m 


BLI 


0.0083444 


0.4 


0.0517 ± 0.0009 


HD149bp41 


E 


BLI 


0.0083444 


0.4 


0.0517 ± 0.0010 


HD149bp42 


n 


BLI 


0.0083565 


0.0 


0.0503 ± 0.0008 


HD149bp42 


i 


BLI 


0.0083564 


4.0 


0.0513 ± 0.0009 


HD149bp42 




BLI 


0.0083562 


5.6 


0.0502 ± 0.0016 


HD149bp43 




BLI 


0.0083681 


0.0 


0.0536 ± 0.0008 


HD149bp43 




BLI 


0.0083682 


6.5 


0.0525 ± 0.0010 


HD149bp43 




BLI 


0.0083685 


7.0 


0.0516 ± 0.0011 


HD149bp43 


m 


BLI 


0.0083682 


7.2 


0.0527 ± 0.0010 



4.2. Joint Fit 

We perform two joint-model fits, each requiring less than 
2 X 10^ iterations to estimate uncertainties. The first consid- 
ers only the three transits analyzed here while the second also 
consid ers the more precise NICMOS data from ICarter et al.l 
(^2009^ by placing priors on / and a/R^. Both fits in Table |5] 
are consistent with previous resu lts from Knutson et al. (201ii 
Rp/R^ = 0.0522 ± 0.0008) and ICarter et al.1 (i200'9l Rp/R^ = 
0.0519 ± 0.0008) and have improved estimates of the ra- 
dius ratio. The uncertainties in the duration and ingress/egress 
times for the independent fit are significantly larger than those 
from the fit with .Carter et al.. (.2009.) priors. 



TABLE 5 
Joint Transit Model Fits 



Paraineters 


Independent Fit 


Carter et al. (2009) Priors 


Rp/R. 


0.0514 ±0.0006 


0.0518 ±0.0006 




87.2+1-^ 


84.6 ± 0.5 


a/R* 




5.98 ±0.17 


Impact Parameter 




0.57 ± 0.04 


Transit depth [%] 


0.264 ± 0.006 


0.268 ± 0.006 


Duration [ti-t4, hr] 


a 9a+0.04 
-^■^ -0.02 


3.286 ±0.019 


Ingress/Egress [hr] 


o.i78!i]:!]« 


0.234 ±0.012 


HD149bp4I 






Midpoint (MJDutc)'' 


4327.3719 ± 0.0005 


4327.3720 ± 0.0005 


Midpoint (MJDxdb)'' 


4327.3726 ± 0.0005 


4327.3727 ± 0.0005 


O-C (minutes)'' 


-1.0 ±0.7 


-0.9 ± 0.8 


SDNR 


0.0083440 


0.0083440 


HDI49bp42 






Midpoint (MIDutc)*" 


4356.1316 ±0.0005 


4356.1316 ±0.0005 


Midpoint (MIDjub)'' 


4356.1323 ± 0.0005 


4356.1323 ±0.0005 


O-C (minutes)'^ 


0.2 ± 0.7 


0.0 ± 0.8 


SDNR 


0.0083550 


0.0083556 


HD149bp43 






Midpoint (MIDuj^)'' 


4597.7070 ± 0.0004 


4597.7068 ± 0.0005 


Midpoint (MIDjub)'' 


4597.7077 ± 0.0004 


4597.7075 ± 0.0005 


O-C (minutes)'^ 


0.8 ± 0.7 


0.6 ± 0.6 


SDNR 


0.0083690 


0.0083691 



"We place priors on (' and a/R^ using values from lCarter et alj 120091) . 

''MJD = BJD - 2,450,000. 

'^Computed using the period and ephemeris from lKnutson et al] 120091 p = 
2.8758925 ± 0.0000023 days, = 2454597.70645 ± 0.00018 BJDutc)- 



5. SECONDARY-ECLIPSE FITS AND RESULTS 

There were 11 secondary-eclipse observations. We con- 
sidered whether the eclipse duration is consistent with the 
more-precise transit duration in these fits, which is likely if 
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the orbit is nearly circular. In a joint fit of all secondary- 
eclipse events, the strong signal from HD149bsll dominates 
the shared eclipse duration. The best-fit eclipse duration was 
4.5 ±3.3 minutes l onger than, but still consistent with, the 
transit duration from lCarter et alJ (12009), and the mid-eclipse 
phases were in all but one case within l.Scr of 0.5, together 
indicating circularity. Since the transit and eclipse durations 
are consistent, we apply priors to the eclipse duration and 
ingress/egress times using the values given in Table|5] Unless 
otherwise stated, we estimated the background flux using an 
annulus from 7 to 15 pixels that was centered on the star We 
present the RMS model residuals in Figure |9] the best-fit light 
curves in Figure [TOl and the best-fit parameters in Appendix 
IB] The electronic supplement contains light-curve files and 
supplementary figures. Below, we discuss each observation 
in detail. 
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Fig. 9. — Same as Figure [7] except for eleven HD 149026b eclipses. The 
shorthanded legend labels correspond to the last three characters in each 
event's label (e.g., sll = HD149bsll). 



5.1. Fit at 3.6 ^im- HD149bsl 1 

For this data set, we find that BLI outperforms NNI down 
to a bin size of 0.015 pixels when we exclude bins with less 
than 4 measurements. Bins with fewer data points are insuffi- 
ciently sampled to compute a reliable mean flux for the knot 
value. The linear ramp (Eq. [TOl i fits best. The posterior pa- 
rameter distributions (histograms) and parameter correlation 
plots (including knot values in the BLISS map) are in Figures 
[TT|-[T3] Similar plots at other wavelengths are included in the 
electronic supplement. 



5.2. Fit at 4.5 /xm - HD149bs21 

Using the strategy described in Section [331 BLI achieves a 
better fit than NNI with bin sizes of 0.012 pixels along x and 
0.006 pixels along y. Additionally, we fit the intrapixel sensi- 
tivity with three different polynomial models ranging between 
second and sixth order After clipping the first 5,000 points 
iq = 5,000), the eclipse depths using various ramp and in- 
trapixel model components are consistent (see Table|6]i; how- 
ever, the SDNR clearly favors BLI and the BIC favors Eq. |2L 
to model the systematics. To minimize the convergence time 
in our MCMC chains, we orthogonalize the eclipse depth, sys- 
tem flux, and both ramp parameters (ro and ri). All of the 
model fits exhibit, to various degrees, bimodal distributions 
in the eclipse-midpoint histograms. The lesser peak occurs at 
a phase of 0.497 and is likely a result of the model trying to fit 
the eclipse egress to the points from phases of 0.514 to 0.520, 
which are consistently above the secondary eclipse hy \ -2a 
(see Figure [TOl center panel). 



TABLE 6 

HD149BS21 - Comparing Model Fits 



R{t) 


M(.v,y) 


SDNR 


ABIC 


Eel. Depth 
[%] 


EL 


BLI 


0.0038800 


0.0 


0.034 ± 0.006 


m 


BLI 


0.0038800 


1.4 


0.033 ± 0.007 


a 


BLI 


0.0038800 


10.8 


0.033 ± 0.006 


a 


Quad. Poly. 


0.0038961 




0.032 ± 0.006 


a 


Cubic Poly. 


0.0038954 




0.033 ± 0.006 


a 


Sextic Poly. 


0.0038941 




0.034 ± 0.006 



5.3. Three Fits at 5.8 /iw 

Using the processes described below, we choose the best-fit 
model for each event, then perform a joint fit with a single 
eclipse-depth parameter. 

5.3.1. HD149bs31 

Initially reported by IStevenson et alj (120101) . we find clear 
evidence of pixelation at 5.8 |^m (see Figure|2|. A relatively 
large bin size of 0.04 pixels is appropriate for HD149bs31 
using BLI in combination with Eq. ^ to express the time- 
dependent flux variation. We orthogonalize the system flux 
and both ramp parameters when computing uncertainties. The 
eclipse-midpoint histogram peaks at a phase of 0.502 and has 
a broad uncertainty of 0.005. The best-fit eclipse depth is 
0.044 ± 0.016%. 

5.3.2. HD149bs32 

The pixelation effect in HD149bs32 is weak because stellar 
centers fall predominantly near the middle of an interpolated 
subpixel (i.e., away from the blue peaks in the right panel of 
Figure|5]l. As such, an intrapixel model component is unnec- 
essary and we model its systematics solely by Eq. |2]^. We 
orthogonolize the same parameters as with HD149bs31. The 
eclipse-midpoint histogram shows a strong bimodal distribu- 
tion, with peak phase values of 0.496 and 0.501. The latter is 
favored with an eclipse depth of 0.038 ± 0.017%. This data 
set points to the illegitimacy of fixing one of the ramp pa- 
rameters (see Section l2.5.3l l. The eclipse depth is correlated 
with the ramp parameters, so fixing one of them erroneously 
improves the eclipse depth uncertainty to 0.012%. 
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Fig. 10. — Binned (left) and systematics- corrected (center & right) secondary-eclipse light curves of HD 149026b in five Spitz.er channels. The results are 
normalized to the system flux and shifted vertically for ease of comparison. The colored lines are best-fit models and the error bars are \a uncertainties. The 
shorthanded legend labels con'espond to the last three characters in each event's label (e.g., si 1 = HD149bsll). 
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Fig. 11. — Parameter histograms for HD 149bs 1 1 . We plot every 200* step 
in the MCMC chain to decorrelate parameter values. The BLISS map knots 
are similarly distributed. Additional histograms are part of the electronic 
supplement. 
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Fig. 12. — Parameter coiTelations for HD149bsll. The background color 
depicts the absolute value of the correlation coefficient. The uncertainties 
produced from the MCMC method fully account for correlations between 
free parameters (e.g., eclipse flux ratio and system flux). We plot every 200* 
step in the MCMC chain to decorrelate parameter values. Additional correla- 
tion plots are part of the electronic supplement. 
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Pixel Position in x 

Fig. 13. — Correlation coefficients between eclipse depth and computed 
BLISS map knots for HD149bsl 1. The presence of relatively strong congela- 
tion regions (in red) indicates that computing the BLISS map at each step of 
an MCMC routine is necessary to assess the uncertainty on the eclipse depth 
correctly, as opposed to fixing the map as is done by BIO. In this case, fix- 
ing the BLISS map to its post-minimizer values leads to an erroneous 13% 
decrease in the eclipse-depth error estimate. 

5.3.3. HD149bs33 

For the same reasons as with HD149bs32, no intrapixel 
model component is needed with HD149bs33. We model 
the declining flux using Eq. [J]., orthogonolize the same three 
parameters as above, and find an eclipse depth of 0.047 ± 
0.016%. The histogram of eclipse midpoints is clearly bi- 
modal but favors the best-fit value of 0.500:!:" Qog- 

5.3.4. 5.8 Joint Fit 

To improve S/N on the eclipse depth, we share this param- 
eter in a joint fit of all three data sets. We retain individ- 
ual eclipse-midpoint times for the subsequent orbital analy- 
sis. The MCMC chain converged after 6 x 10^ iterations. The 
combined light curve in Figure[T4]illustrates the improvement 
when compared to the three 5.8 [a.m light curves in Figure [TOl 
The best simultaneous fit favors an eclipse depth of 0.044 ± 
0.010%. 

5.4. Four Fits at 8.0 fim 

Similarly to the fits at 5.8 \xm, we fit models to each data set 
individually then use the best models in an 1.1 x 10^ iteration 
joint fit that shares a common eclipse depth. 

5.4.1. HD149bs41 

The significant improvements in our pipeline since the orig- 
inal analysis by H07 warrant a new analysis of HD149bs41. 
We follow all of our current techniques described in Section|2] 
and test all of the listed ramp model components. As with H07 
and K09, we find that Eq. |2]^ best describes the overall ramp. 
The smaller ramps associated with each telescope movement 
are best described by Eq. [13] according to BIC; however we 
also present H07's 12-point spline for comparison (see Table 
|7]l. Each model employs a constant-flux offset at each of the 
nine nod positions, of which eight are free parameters as de- 
scribed in Section|2] 

Due to the nodding motion with this particular data set, BLI 
and NNl are inappropriate models to use. We can see in Fig- 
ure[T5] which illustrates one of the nine nod positions, that the 




0.9985 
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Fig. 14. — Combined and binned eclipse light curves at 5.8 and 8.0 |^m. 
No te th e improved S/N achieved with combined modeling, compared to Fig- 
ure[T0] The best joint-fit HDI49bs32 and HD149bs41 models are plotted for 
comparison. 



pixel position is slightly different for each of the twelve visits 
to this position. This behavior introduces a strong time de- 
pendence in the position sensitivity correction that cannot be 
disentangled. Our best attempt to correct for the position sen- 
sitivity uses a linear correction in two dimensions for each of 
the nine nod positions (9 x Linear). Table|7]compares the four 
best model combinations. Compared to K07, our flux offsets 
are multiplicative rather than additive (see Eq.[T]i and our final 
model does not fix either of the ramp parameters. As a result, 
our eclipse-depth uncertainty is larger (0.049 ± 0.016%). 




17.30 17.35 17.40 17.45 17.50 17.55 

Pixel Position in x 

Fig. 15. — Pointing histogram for one of nine nod positions of HD149bs41. 
The small, 0.01-pixel bin size clearly shows that the positions of the 12 visits 
have very little overlap, resulting in a time-dependent position sensitivity and 
making the data impossible to model accurately using a BLISS map. The 
small footprint size demonstrates the difficulty of making a definitive IRAC 
intrapixel map using all the stellar staring data in each channel. The horizon- 
tal and vertical black lines represent pixel boundaries. 
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TABLE 7 

HD149BS41 - Comparing Model Fits 



R(t) 


M(x,y) 


Visit Sensitivity 


SDNR 


ABIC 


Eel. Depth 

[%] 


E- 




Eq.\n\ 


0.0084575 





0.049 


EL 




12-pt Spline 


0.0084559 


59 


0.049 


EL 


9 X Linear 


Eq.[T3] 


0.0084510 


126 


0.050 


EL 


9 X Linear 


12-pt Spline 


0.0084494 


184 


0.050 



5.4.2. HD149bs42 

We use a 0.04-pixel bin size and only consider bins with 
at least eight points to ignore an outlier near subpixel loca- 
tion (14.36, 15.24). Table [8] contains ABIC values and best- 
fit eclipse depths from our least-squares minimizer for three 
different values of the clipping parameter: q=Q, 5,000, and 
10,000. This parameter ignores the given number of data 
points fror n the beginning of the observation and is a common 
procedure (iKnutson et al.ll20TT ') when trying to find the best- 
fitting ramp. The table indicates that all but one of the eclipse 
depths are consistent for q = 10,000 and that Eqs.|4]_ and|5l_ 
are consistent at all three q values. Our MCMC routine finds 
strong non-linear correlations in all ramp model components 
except Eq. [TT] which exhibits linear correlations that are eas- 
ily handled by orthogonalizing the system flux and both ramp 
parameters. We use Eq. [TTIwith q = 10,000 in the joint model 
fit. The eclipse depth for the competing solution (Eqn. |2]_) 
differs by less than \a. 



TABLE 8 

HD149BS42 - Comparing Model Fits 



q: 

R(t) 




Eel. Depth 

[%] 




ABIC 


5000 
Eel. Depth 
[%] 


5000 
ABIC 


10000 
Eel. Depth 

[%] 


10000 
ABIC 




0.149 


59 


0.111 


1 


0.068 







0.030 


4 


0.056 


2 


0.068 


11 


w 


0.065 


7 


0.065 


9 


0.063 


19 


S- 


0.065 


8 


0.062 


9 


0.061 


20 


m 


0.132 


47 


0.099 


6 


0.067 


10 


m 


0.096 


37 


0.086 


15 


0.065 


20 




0.087 





0.071 





0.068 


10 


m 


0.062 


16 


0.052 


17 


0.049 


28 


m 


0.197 


205 


0.128 


31 


0.064 






5.4.3. HD149bs43 

We choose a relatively large bin size of 0.03 pixels because 
the intrapixel effect is minimal (position dependence is flat) 
and we do not want to overfit the edges of position space 
where there are few data points. Smaller bin sizes result in 
similar eclipse depths. Without BLl, we find an eclipse depth 
of 0.040 ± 0.008%, nearly identical to K09; however, this fit 
does not have the lowest BIC value (ABIC = 1 12). There is a 
small but clear position dependence that, once accounted for, 
results in a marginally deeper eclipse of 0.044 ± 0.008%. We 
confirm that the eclipse midpoint is noticeably earlier (by Aa) 
than the expected phase value of 0.5, but do not claim a detec- 
tion of eclipse timing variation because HD149bs21 measures 
the preceding eclipse with a much stronger S/N and is consis- 
tent with a circular orbit. Fixing the phase of mid-eclipse to 
0.5 results in a marginally shallower eclipse depth (0.039 ± 
0.008%) and a larger SDNR value. 



5.4.4. HD149bs44 

Relative to the other 8.0 |J.m events, HD149bs44 exhibits 
significantly larger SDNR and uncertainty scaling factor val- 
ues (see Table[T6lin AppendixlBj. The uncertainty scaling fac- 
tor renormalizes the error bars such that = 1, so a smaller 
scaling factor indicates a better fit. For this noisy data set, the 
models achieve relatively poor fits compared to other 8.0 |J.m 
data sets. 

The ramp models that produce the most consistent eclipse 
depths in Table |9] are Eqs. gL and E] for ^ = 5,000, 7,500, 
and 10,000. Smaller values of q produce inconsistent eclipse 
depths for all ramp models; larger q values do not provide 
sufficient out-of-eclipse baselines. For q = 7,500, most of the 
ramps find eclipse depths that are in agreement with the con- 
sistent values given by Eqs.HL and[TT] Of these ramps, Eqs. 
|2L and [TT] share the lowest BIC value; however, the quadratic 
parameter in Eq.[TT]correlates strongly with the eclipse depth, 
resulting in a larger uncertainty (0.017 vs. 0.013), so we select 
Eq.|2]_ with q = 7,500 for our final joint model. This data set 
also applies a BLISS map with 0.025-pixel bin sizes and at 
least ten points per bin. 



TABLE 9 

HD 149BS44 - Comparing Model Fits 



q: 5000 


5000 


7500 


7500 


10000 


10000 


R(t) Eel. Depth ABIC Eel. Depth ABIC Eel. Depth ABIC 


[%] 




[%] 




[%] 




El 0.072 





0.066 





0.085 





[E 0.075 


11 


0.069 


11 


0.085 


11 


|4]_ 0.073 


23 


0.068 


22 


0.069 


24 


m 0.072 


22 


0.066 


22 


0.085 


21 


|6] 0.074 


11 


0.069 


11 


0.082 


11 


|7] 0.089 


17 


0.096 


10 


0.086 


22 


m 0.073 


33 


0.067 


22 


0.081 


32 


[ni 0.073 


1 


0.069 





0.069 


3 



5.5. Two Fits at 16 fim 

Each event consists of 1050 exposures, divided serially into 
three 350-image sequences at non-overlapping stellar centers 
on the detector The second sequence ip=l) is completely 
within the secondary eclipse and, because of the free position 
offset parameter, has no impact on the eclipse depth (but still 
helps to model the systematics). Unfortunately, this means 
that the eclipse depth is completely determined by the small 
number of points after ingress in the first sequence and before 
egress in the last sequence. 

To avoid residual effects from potentially bright previous 
targets, neither of the two 16 \xm observations was positioned 
at the center of the array. As a result, the outer radius of the 
sky annulus is limited to 1 1 and 12 pixels for HD149bs51 and 
HD149bs52, respectively, to avoid the flux falloff at the edges 
of the blue peak- up array. We apply both aperture and opti- 
mal photometry (lHornelll986l; iDeming et al]|2005h ; flie for- 
mer produces cleaner results for both data sets. 

Figure [16] displays the best-fit eclipse depths and corre- 
sponding SDNR values versus aperture size for four compet- 
ing models in each data set. The first applies no intrapixel 
correction, the second uses three linear components (one at 
each position), the third uses BLl, and the final model uses 
NNl. The first two models also apply a position offset to ac- 
count for the differences in pixel sensitivity. This offset is 
redundant for BLl and NNl. 
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Although these models find similar (< Icr) eclipse depths at 
the best aperture size of 2.0 pixels, this is not the case for other 
aperture sizes. Models using BLl or NNl produce the most 
consistent eclipse depths and result in lower SDNR values 
(see Tables fTOlandfTTTl. We test a range of bin sizes for BLISS 
mapping but find that NNl consistently outperforms BLl. We 
take this as evidence that there are insufficient data points in 
most bins to model the weak position dependence accurately. 
As such, we consider only the last two models and select 
no intrapixel model for HD149bs51 and the 3 x linear model 
for HD149bs52 to perform the joint fit. For improved con- 
vergence in our MCMC chains, we orthogonalize the system 
fluxes and both ramp terms in Eqs.[TT]and|2]_ for HD149bs51 
and HD149bs52, respectively. Due to the weak eclipse signal, 
the midpoint is fixed to a phase of 0.5 for the individual fits; 
the joint fit achieves a sufficient S/N to share the eclipse mid- 
point as a free parameter. The final model fit uses 1.6 x 10^ 
iterations and is in Table [T6]of Appendix IE] 



HD149BS51 



TABLE 10 

■ Comparing Model Fits 



R(t) 


M(x,y) 


SDNR 


ABIC 


Eel. Depth [%] 


\m 


NNl 


0.002890 





0.059 ± 0.038 


m 


BLl 


0.002924 


21 


0.061 ± 0.035 


m 




0.003119 


169 


0.068 ± 0.038 


m 


3 X Linear 


0.003066 


175 


0.060 ± 0.037 






TABLE 1 1 






HD149BS52 - Comparing Model Fits 


R(t) 


M(x,y) 


SDNR 


ABIC 


Eel. Depth [%] 




NNl 


0.003166 





0.081 ± 0.036 




BLl 


0.003211 


23 


0.068 ± 0.035 


EL 


3 X Linear 


0.003378 


175 


0.074 ± 0.038 


EL 




0.003574 


243 


0.066 ± 0.034 



We find that the eclipse depths for both data sets vary with 
the choice of eclipse midpoint. Using the midpoint from our 
joint model (0.5015, see Table [TSTi, we find that the best-fit 
eclipse depth differs by up to la from the values listed in Ta- 
bles [TO] and (TT] where the midpoint is fixed to 0.5. We also 
test joint fits without the data at p = 1 and find a comparable 
joint-fit eclipse depth of 0.099 ± 0.035%. The small change 
is likely the result of weak correlations with the ramp param- 
eters. We include all positions in the final fit as this results in 
a lower SDNR and a small improvement in the eclipse depth 
uncertainty. 

6. ATMOSPHERE 

In order to understand the atmosphere of the planet bet- 
ter, we compare our measured flux ratios to those generated 
from a model atmosphere. We simulate t he atmosphere of HD 
149026b using the model presented by iFortnev et al.l (I2005L 
|2006, 2008). See IFortnev et al.l (l2008b for a description of 
the heritage of the model, which includes solar system plan- 
ets and brown dwarfs in addition to exoplanets. The chemical 
mixing ratios used assum e chemical equilibrium, following 
iLodders & FeglevI (120021) . at both solar metallicity ("Ix sO' 



lar") a s well as 30 x solar, using the abunda nces of ILodders 



(2003"). The opacity database is described bv lFreedman et al 
(12008.) . with an update to include CO2 opacity. 
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Fig. 16. — Best-fit eclipse depths and eon'esponding SDNR values for 
HD149bs51 and HD149bs52. Both are plotted as a function of photome- 
try aperture size for the HD149bs51 (upper panel) and HD149bs52 (lower 
panel) data sets. The labels refer to the type of intrapixel model component 
used. Here, 'None' uses no model, '3*Ln' uses a linear function in .v and y 
at each of the three positions, and 'BLF and 'NNF both use our new BLISS 
mapping technique with 0.02-pixel bins. In all but one case, the best aperture 
size is 2.0 pixels, according to SDNR. Due to the weak eclipse signal, the 
phase of mid-eclipse is fixed to 0.5. A typical Icr eclipse-depth uncertainty is 
0.037%. 



We have generated chemistry/opacity grids with and with- 
out the opacity of gaseous TiO/VO. These gases, which are 
strong absorbers of optical flux, may be responsible for the 
temperature inversions diagnosed in the atmosphere of some 
planets (e.g., Hubeny et al. 2003; Fortney et al. 2006, 2008), 
but see Spiegel et al., (,2009) for important caveats. The mid- 
infrared flux ratios observed for hot Jupiters are quite di- 
verse, but high flux ratios (and corresponding large brightness 
temperatures) in the mid-infrared, together with small 3.6-to- 
4.5 |.im rati os, have been found in models wi th temperature 
inversions (IFortnev et al.l I2006| iB urrows et al. 20 07., 20081 
iMadhusudh an & Seager 2009; Spiegel & Burrows i2010l) . and 
thes e models have had some success in comparisons with data 
(see lSeageF & Deming 2010, for a review). For HD 149026b, 
we find a 3.6-to-4.5 [a.m ratio > 1, similar to HD 189733b, 
which indicates no temperature inversion. 

We compare the measured flux ratios to three different mod- 
els in Figure [it] All assume redistribution of absorbed stellar 
flu x over the dayside of the planet only (/ = 1/2, as described 
by IFortnev & Marlevll2007 '). We show a 1 x solar model with 
TiO/VO opacity (which yields an inversion) and the same 
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model with TiO/VO opacity neglected. We also plot a sim- 
ilar no-inversion model at 30 x solar metallicity. This last 
model leads to mixing ratios ^-^30 and ^900 times larger for 
CO and CO2, respectively, compared to the solar metallicity 
case ( Zahnle et al.li2d09h . The dramatic increases in CO and 
CO2 lead to the most notable spectral difference, the much 
deeper absorption due to the overlapping 4.5 ^.m band of CO 
and 4.2 \im. band of CO2. Such a high metallicity for this 
atmospher e may well be rea listic, as Saturn is ~ lOx solar 
in carbon ( Flasar et al.ll2005h . while Uranus and Neptune are 
~ 30-60 X enriched. 

HD149026b Day 
0.001 2 r"' ' ' ' \ ' ' ' ' I ' ' ' ' I ^ 




Fig. 17. — Atmospheric models of the dayside of HD 149026b. The red 
line depicts a model with a temperature inversion, / = 0.50, solar metallicity, 
and an effective day-side temperature, Tgff, of 2067 K. The blue model has 
no temperature inversion, / = 0.50, solar metallicity, and T^ff = 2056 K. The 
green model is similar to the blue model but with enhanced metallicity (30 X 
solar) and T^ff = 208 1 K. Model band-averaged ratios are shown as squares 
while the data points are orange diamonds. Models that lack temperature 
inversions and include a high atmospheric metallicity best match the observed 
data. 

Models with a temperature inversion similar to the red line 
shown in Figure[T7]are clearly disfavored, given the 3.6-to-4.5 
|j.m ratio > 1 and the large flux ratios at redder wavelengths. 
Our best fit for the three models is the 30 x enhanced, no- 
inversion model. Only the 3.6 |xm point is outside our la error 
bar Even with the strong CO/CO2 absorption, we still cannot 
quite match the observed 3.6-to-4.5 \un ratio. At face value, 
this would imply even larger CO and/or CO2 mixing ratios; 
however, this may not necessarily be the case, as some other 
modelers (e.g. Burrows et al. 2008) generally show a deeper 
absorption feature at 4.5 |J.m than we obtain with our models. 
This is likely due to differences in the temperature gradient as 
a function of height in the different models, as this helps to 
control the depth of absorption features. A steeper gradient 
leads to deeper absorption features. 

Less-efficient redistribution of absorbed flux leads to a hot- 
ter dayside, and still would yield a satisfactory fit to the ob- 
servations, albeit near the top of the Icr error bars. Given the 
8.0 |a.m phase curve of K09, which showed modest day/night 
phase variation, such a hot dayside (which would leave little 
energy for the night side) is not favored. We recommend that 
a coupled 3D dynamics/radiative transfer model be run for the 
system, to understand if the imphed day- and nigh t- side tem- 
peratures can be matched. IShowman et al.l (l2009t) had good 
success in matching the phase curve and dayside photometry 
of HD 189733b. These models tend to show better day-night 
homogenization of temperature contrasts than one would as- 



sume from the fact that our best-fit ID model assumes all ab- 
sorbed flux is re-radiated on the planet's day side. Addition- 
ally, n ear-infrared fluxes from the JHK-band (e.g., ICroll et alj 
120 101) , where this planet is brightest, would help to understand 
the dayside luminosity; however, given the small planet-to- 
star radius ratio, this may have to wait for the James Webb 
S pace Telescope . 

iKnutson et alj (12010 1) suggest that the absence of a tempera- 
ture inversion within an exoplanet atmosphere correlates with 
higher levels of chromospheric activity from the host star. The 
lack of a temperature inversion in HD 149026b does not agree 
with HD 149026's relatively low activit y level, but t his may 
be due to the exoplanet's high density (.Knutson et aLllTolQ) . 

The pressure-temperature profiles for the three atmospheric 
models are shown in Figure [18] Also plotted are the contri- 
bution functions (e.g.. Chamberlain & Hunten 1987) for ther- 
mal flux in each of the five Spitzer bandpasses. Contribution 
functions trend towards lower pressure with enhanced metal- 
licity for all bandpasses, but move most dramatically at 4.5 
|j.m due to CO and CO2. At constant metallicity, a temper- 
ature inversion tends to smear the contributions over a wider 
pressure range, favoring lower pressures due to the hot upper 
atmosphere, but to a lesser extent compared to models with 
increasing metallicity. 

7. ORBIT 

We have collected a total of eleven individual Spitzer 
secondary-eclipse observations with useful timing of HD 
149026b over a 3.5 year baseline. These times constrain 
ecosw, where e is the eccentricity and cj is the argument of 
periapsis, and can be used to establish eccentricity limits on 
the planet's orbit. We use BJDjQg given in Tables [TSl and [T6l 
and correct for the eclipse-transit light-time (42 seconds). The 
mean eclipse phase, using the K09 ephemeris, is 0.49997 ± 
0.00028, suggesting that ecosw = -0.00003 ± 0.00044. The 
data are consistent with a circular orbit (e < 0.0013). The 
times of secondary eclipse do not show any significant trends 
and do not have a period that differs significantly from the pe- 
riod determined from transit and radial velocity data. Such a 
differ ence would indicate apsidal motion or other secula r ef- 
fects ( Gimenez & Bastero|[l995l;lHevl & Gladmanll2007 ^. An 
MCMC ephemeris fit to our secondary-eclipse times gives a 
period of 2.875884 ± 0.000006 days, and a fit of all the avail- 
able transit times gives a period of 2.8758922 ± 0.0000015 
days. The difference between the two periods is not signifi- 
cant ( 1 Act). If the measured period difference is due to apsidal 
motion, then it would indicate that cOesina; = (9 ± 7) x 10"^ 
°/day ( Gimenez & Bastero 1995), where lu is the rate of ap- 
sidal precession. Further secondary-eclipse observations will 
refine the secondary-eclipse period. 

We use our primary-transit and secondary-eclipse data to 
perform a fit, as described by Campo et al. (201 1), that also 
incorpo rates other available transit data ( C arter et al.l 120091 ; 
iWinn et al. .2008; .Charbonneau et al. .200 6) and radial veloc- 
itv data dSato et al.ll2005l; iButler et aljl200 6). When we fit an 
eccentric orbit to the available data (Table [T2]l, we determine 
that e = 0.154 ±0.016. Although this is a 1 Oct eccentricity, it is 
almost completely dominated by the e sin lj component, lead- 
ing us to believe that the eccentricity may be an overestimate 
(Laughlin et al. 2005). This is possible when the peaks of the 
radial velocity curve, where the waveform is most sensitive to 
changes in esinw, are undersampled. The eccentricity affects 
the symmetry of the RV curve, so when both peaks are not 
well sampled, the best-fit solution may misrepresent the ac- 
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Fig. 18. — Contribution functions a nd a tmospheric pressure-temperature profiles. Left 3 panels: The contribution functions in the five Spitzer bandpasses are 
for the three models shown in Figure [TT] Emission generally comes from higher in the atmosphere for the metal-enriched model (left) and, to a lesser extent, 
the model th at fe atures a temperature inversion (right). Right-most panel: The atmospheric pressure-temperature profiles are for these same models, colored to 
match Figure [TtI The 30 X solar no-inversion model is everywhere warmer than the 1 X solar no-inversion model, but they have very similar T^ff values, since the 
emission from the 30 X model comes from much higher in the atmosphere. 



tual eccentricity of the planet's orbit. Indeed, 16 of the 23 us- 
able RV data points were taken at a transit phase greater than 
0.5 and there are no data points between phase values of 0.1 
and 0.3. The dearth of RV measurements near 0.25 signifies 
that only one of the two peaks is adequately constrained. To 
best refine the value of esincj, we require additional RV mea- 
surements between phases and 0.5, particularly near 0.25. 

In a comparison fit assuming a circular orbit (see Table [13), 
where the RV curve is perfectly sinusoidal and symmetric, 
BIC is worse than for the eccentric fit. Despite this, the un- 
dersampled radial velocity data and the high degree of con- 
sistency of the eclipse phases with 0.5 make it unlikely that 
the orbit of this planet has an eccentricity greater than the 
maximum value of |ecosw|. A near-perfect alignment of the 
system's semi-major axis with our line-of-sight {u ^ 90°) 
would be necessary, but the agreement between the transit and 
eclipse durations (see Section |5]l argues against this scenario. 
Acknowledging that our secondary-eclipse timing measure- 
ments yield little information about esinuj, we present both 
solutions without judgment. Although an eccentric orbit is 
unlikely, it cannot be ruled out with the data currently avail- 
able. 

8. CONCLUSIONS 

Over 3.5 years, Spitzer observed three primary transits 
and eleven secondary eclipses of HD 149026b in five in- 
frared wavelengths. We utilize multiple observations for 
channels with the weakest eclipse depths to improve S/N es- 
timates and better constrain the dayside atmospheric com- 



TABLE 12 
Eccentric Orbital Model 



Parameter 


Value 


esinw 


0.154± 0.016 


ecosoj 


-0.00037 ± 0.00044 


e 


0.154 ±0.016 




90. 14 ±0.16 


P (days) 


2.8758919 ±0.0000014 


To (MJDtdb)" 


4597.70716 ±0.00016 


(ms'') 


47.4± 1.1 


7 (ms"') 


-4.3 ±0.6 


BIC 


123 



'MJDxDB = BJDxDB - 2,450,000 



TABLE 13 
Circular Orbital Model 



Parameter Value 
P (days) 2.8758916 ±0.0000014 
To (MJDxdb)" 4597.70713 ± 0.00016 
K{m?,-^) 42.6 ±0.9 
7(ms-l) -1.6 ±0.6 
BIC 179 

"MJDxDB = BJDxDB " 2,450,000 

position. The addition o f a third transit event at 8.0 ixm 
confirms previous results dNutzman et al.l l2009t ICarter et"an 
120091: iKnutson et al.ll20()9l) and offers an improved constraint 
on the planet-to-star radius ratio {Rp/R^). A new ecHpse anal- 
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ysis of HD149bs41 confirms the findings froml Knutson et aLl 
( 120091) . namely that an eclipse depth of ^ 0.05% fits best at 
8.0 i^m. However, we find a larger uncertainty due to corre- 
lations between the eclipse depth and ramp parameters that 
were not fully explored because one of the ramp parameters 
was previously fixed. 

The atmosphere is explained well by a ID chemical- 
equilibrium model. A temperature inversion is no longer 
favored when fitting the observed planet-to-star flux ratios. 
The best-fit model includes large amounts of CO and CO2, 
moderate heat redistribution (/ = 0.5), and strongly enhanced 
metallicity (30x solar). Using the times from our secondary- 
eclipse observations, we find no deviations from a circular or- 
bit at the 1 a level. However, given the available RV data, we 
cannot completely rule out an eccentric orbit with an unlikely 
orbital alignment. 

We present a new technique, called BiLinearly-lnterpolated 
Subpixel Sensitivity (BLISS) mapping, to model Spitzer's 
position-dependent systematics (intrapixel variabiUty and pix- 



elation). In all cases tested to date, BLISS mapping outper- 
forms previous methods in both speed and goodness of fit. 
We also apply an orthogonalization technique for linearly- 
correlated parameters that accelerates the convergence of 
Markov chains that employ the Metropolis random walk sam- 
pler. 
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A. MODEL COMPARISON WITH COMPLEX MODELS 



In Section|3] we show that BLISS m apping impro ves the SDNR of models that include an intrapixel sensitivity term. Elsewhere 
in this paper and in our earlier work (ICampo et a l. 201 1) we use the Bayesian Information Criterion (BIC) to compare models 
of other systematic error components (e.g., rival parametric ramp models), but we did not use the BIC to compare intrapixel 
sensitivity models. In this Appendix we discuss the approximations underlying the BIC to elucidate when it is useful, and in 
particular, to explain why we do not use it (or similar statistical criteria) for comparison of intrapixel sensitivity models. 

The BIC provides an asymptotic approximation to quantities that may be used for Bayesian quantification of model uncertainty. 
The most simple use of the BIC is to approximate Bayesian explanatory model selection. A basic distinction among model 
selection criteria is between explanatory criteria (that seek the model that best describes the processes that produced the available 
data) and predictive criteria (that seek the model that will make the most accurate predictions of future data based on the available 
data). For example, the well-known Akaike Information Criterion (AIC = x^ + 2fe) is a predictive criterion. Although different 
from the BIC in rationale, its derivation invokes similar asymptotic approximations to those we describe here for the BIC, and 
we consider it to be similarly hampered for comparing large models. To illuminate the nature of the approximations underlying 
the BIC, we sketch its derivation here. 

Consider a set of models {M,} for observed data D. Let 6*, denote the parameters of model /, with A:, dimensions. The 
(Bayesian) posterior probability for model / is proportional to the product of a prior probability for the model, and the model's 
marginal likelihood. 



where 7r,(6',) = /?(6',|M,) is the prior probability density function (PDF) for the model's parameters, and £,(6',) = /?(D|0,,M,) is 
the likelihood function, the probability for the data presuming the model holds with parameters 6', (the likelihood function is 
proportional to exp[-x^(0,)/2] for our models). As its name suggests, the role of the marginal likelihood in quantifying model 
uncertainty is completely analogous to the role of the more familiar likelihood function in quantifying parameter uncertainty 
(within a particular model). Note that 7W, is an average of the likelihood function for model /, not a maximum: in Bayesian 
inference, the weight of the evidence for a model is given by the typical value of the likelihood function for its parameters, not 
the optimum (largest) value. 

For nonlinear models with more than a few dimensions, calculation of the integral in Eq. [17] is not feasible using standard 
quadrature methods. The development of algorithms for accurate calculation of marginal likelihoods is an active research area, 
and existing algorithms are typically problem-specific and computationally expensive (see Clyde et al. 2007 for a recent review 
targeting astronomers). 

The BIC approximates -21nA^ for straightforward models in the limit of voluminous data, i.e., asymptotically (we drop the 
model index, /, when referring to a generic model, to simplify notation). In this limit, the likelihood function C{9) will be strongly 
peaked at the maximum likelihood estimate, 6, and it will be much more strongly concentrated than the prior; Figure [T9l depicts 

the situation. As a first step in approximating the integral in equation[T2l we evaluate the prior at 9 and pull it out of the integral, 

so 



The prior PDF has dimensions of [1/6*], and we can express 7t{9) as the inverse of a local prior scale, A9 (for a normalized 
flat prior spanning a range A9, we have Tr(9) = 1 /A9 exactly). We next approximate the integral of the likelihood function in 
Eq.[T8]as the product of its height (the maximum likelihood value, C(9)) and a characteristic width, S9 (describing the posterior 
uncertainty in 9). This gives 



With suitable definitions of the prior and posterior scales, we can make this equation exact (e.g., for a one-dimensional model 
with a Gaussian likelihood function of standard deviation a, and a flat prior with range A9 ^ a, as long as 9 is not near a 
boundary of the prior range, then setting 59 = a\/27r w 1 .06 times the full width at half maximum, makes the equation accurate). 
The factor multiplying the maximum likelihood is sometimes called the Ockham factor and will be < 1; it quantifies how much 
of the parameter space of the model is wasted, in the sense of including parameter values that are ruled out by the data. Note that 
for dimension A: > 1, these scales are volumes, i.e., products of scales in each dimension. 

Asymptotically, for a simple model of fixed dimension, we expect the uncertainty for each parameter to scale like 1 / \/N for 
sample size A^. So for a model of dimension k, we expect 59 to eventually decrease proportional to A^"*/^. Let Na be the sample 
size where the asymptotic behavior kicks in, and 59a be the typical scale of the uncertainties at that sample size. Then we expect 




(17) 




(18) 




(20) 



In the simple case of estimating Unear parameters for data with additive Gaussian noise of fixed noise variance, we expect 



20 



Stevenson et al. 




AO 



e 







Fig. 19. — Ingredients for the derivation of tlie BIC approximation to tlie marginal likelihood. Curves show the likelihood function (solid blue) and prior PDF 
(dashed green), with characteristic widths 59 and AO. Points show the maximum likelihood parameter estimate, 9, and the values of the likelihood and prior at 9. 



asymptotic behavior right away, so A^a = 1 . Now take the logarithm and group terms according to their dependence on A^: 



The first term may have a nontrivial dependence; the last term is constant with respect to A^. Schwarz (1978) derived a 
more precise expression like this, explicitly calculating the first term in the case of linear models with sampling distributions 
in the exponential family (which includes, e.g., normal, Poisson, and multinomial distributions). In that case the maximum 
log-likelihood term is expected to grow increasingly negative, roughly proportionally to A^. 
The BIC keeps the A^-dependent terms in Eq.|2T|and multiplies by -2; 



It is attributed to Schwarz (and sometimes called the Schwarz criterion), but notably, Schwarz did not drop the A^-independent 
term in his approximation to In Al, although he termed it a "residual" with respect to the A^-dependent terms. 

If the models under consideration are considered equally probable a priori, the most probable model is the one with the largest 
marginal likelihood. BIC -based model selection uses the BIC to approximate the logarithm of the marginal likelihood, choosing 
the model with the smallest value of the BIC. The derivation sketched above provides some insight regarding when we might 
expect this procedure to identify the highest likelihood model. There are two main considerations. 

First, the BIC is an asymptotic criterion. Its accuracy requires sample sizes large enough so that the parameter uncertainties 
are decreasing at the 0(N~'^^^) rate. For complex models with many parameters, this is not simply a matter of sample sizes 
being "large enough." For some models — e.g., nonparametric models — the number of parameters may grow with sample size 
(explicitly or effectively); for others, some parameters may be sensitive to only a subset of the data. For example, the BLISS 
model uses a piecewise linear intrapixel sensitivity map, so a particular coefficient is determined by only a subset of the image 
pixels. In these cases, a BIC-like criterion may be valid, but with the A:lnA^ term replaced with a term that more accurately 
describes the asymptotic behavior of parameter uncertainties. Determining the form of such a term can be subtle (see Kass and 
Raftery 1995, § 4.2). In these settings, it may take very large sample sizes to reach asymptotic behavior. 

Second, the BIC drops a constant (in A^) term from the logarithm of the marginal likelihood — Schwarz's "residual." That 
is, it drops a multiplicative factor from the estimated marginal likelihood. This factor depends on the prior volume, A6. For 
models with many similar degrees of freedom, like the various intrapixel sensitivity models, the prior volume is the product of 
the ranges of many variables. It can vary sensitively with the choice of a priori scale per parameter, and if the competing models 
have different types of parameters, the omitted residual terms may be very different from one model to another The difference 
between the residual terms can be large when the models are large. As a result, the change in the BIC between two large models 
cannot be relied upon for identifying the model with the larger marginal likelihood. 

Kass and Wasserman (1995; KW95) have examined the role of the residual term in the asymptotic approximation of the log 
marginal UkeUhood, arguing that for some problems there may be a reasonable argument for it to be negUgible. Note that the last 




(21) 



BIC = -2 In C(9) + klnN= + klnN. 



(22) 
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term in Eq.|2T]will vanish if 



(23) 



When the asymptotic sample size Na = I (e.g., for linear models, like estimating the mean of a normal distribution with known 
variance), this requirement corresponds to having the prior range equal to the width of the likelihood for a single-sample dataset. 

For Na > the Nc/~ factor scales the S9a likelihood volume to what the single-sample volume would be if the model were 
asymptotic starting with N = I. Thus KW95 dubbed a prior satisfying Eq. |23]a unit information prior, i.e., a prior that is as 
informative as a single datumQ To the extent that one could consider the uncertainty scale associated with a single measurement 
to reflect prior uncertainty, the BIC may be an accurate approximation to InA^. 

However, even when a unit information prior scale appears reasonable, for large models, even a small variation from this scale 
for the prior range of each parameter could produce a large net residual term. We thus do not consider the unit information prior 
argument to provide a sound justification for using the BIC to compare large models. 

For these reasons, in our work we limit use of the BIC to comparing small models, or large models that are nested, so that rival 
models share the vast majority of parameters. For example, we rely on the BIC to compare different ramp models that share a 
common BLISS map model, but we do not consider the BIC to be valid for comparing, say, a model using a BLISS map to a 
model using BlO's intrapixel variability correction, or to polynomial intrapixel models. 

The residual issue, in part, reflects an inherent weakness of marginal-likelihood based model comparisons with large models. 
Small changes in the per-parameter prior scale for such models can lead to large changes in marginal likelihoods, even with an 
accurate numerical calculation of the marginal likelihoods. This motivates developing a way to allow the scale to adapt to the 
data. In a Bayesian framework, this can be accomplished using a hierarchical model to implement regularization. One considers 
the prior range to be an adjustable parameter itself that is learned from the data. This can be done in a manner that essentially lets 
the data determine the effective number of free parameters (the total number of knots) required for the intrapixel map. We would 
then compare our best fit using BLISS mapping to those obtained using other intrapixel models and select the appropriate model 
component. Such calculations are beyond the scope of this paper, but would be a productive avenue for future research. 



^ The KW95 derivation is of course more careful than tliat sketched here. 
They account for correlations between parameters, replacing Eg. 1231 with a 



relationship between Hessian matrices of the prior and likelihood. 
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B. BEST-FIT PARAMETERS 



TABLE 14 

Best-Fit Joint Transit Light-Curve Parameters 



Parameter 


HD149bp41 


HD149bp42 


HD149bp43 


Wavelength (|xm) 


8.0 


8.0 


8.0 


Array Position (x, pix) 


14.93 


14.96 


15.14 


Amy Position (y, pix) 


15.14 


14.56 


14.47 


Position Consistency' (S^, pix) 


0.013 


0.011 


0.011 


Position Consistency' (<5y, pix) 


0.012 


0.013 


0.013 


Aperture Size (pix) 


3.5 


3.5 


3.5 


Sky Annulus Inner Radius (pix) 


7.0 


7.0 


7.0 


Sky Annulus Outer Radius (pix) 


15.0 


15.0 


15.0 


System Flux Fs (\ily) 


1 17229 ± 14 


117460 ±60 


1 17356 ± 13 


Transit Midpoint^ (MJDutc) 


4327.3720 ± 0.0005 


4356.1315 ±0.0005 


4597.7067 ±0.0005 


Transit Midpoint^ (MJDxdb) 


4327.3727 ±0.0005 


4356.1323 ±0.0005 


4597.7075 ± 0.0005 




0.0518 ±0.0006 


0.0518 ±0.0006 


0.0518 ±0.0006 


cos i 


0.095 ± 0.009 


0.095 ± 0.009 


0.095 ± 0.009 


o/R* 


5.98 ±0.17 


5.98 ±0.17 


5.98 ±0.17 


Ramp Equation (R{1)) 


S- 


1 


13- 


Ramp, ro 


24 ±4 





18±2 


Ramp, ri 


-0.6 ±0.7 





4.1 ± 1.4 


Ramp, r2 


0.01 10 ±0.0015 








Ramp, rg 





0.0009 ± 0.0005 





Ramp, rj 





-0.00048 ±0.00011 





Ramp, to 





-0.0248 ±0.0011 





BLISS Map(M(A-,.v)) 


Yes 


Yes 


Yes 


Min. Number of Points Per Bin 


4 


4 


4 


Total Frames 


67008 


54080 


70000 


Rejected Frames (%) 


0.714840 


0.377219 


0.504286 


Frames Used^ 


61520 


53865 


68646 


Free Parameters 


8 


5 


4 


AlC Value 


184044 


184044 


184044 


BIC Value 


184216 


184216 


184216 


SDNR 


0.0083440 


0.0083556 


0.0083691 


Uncertainty Scaling Factor 


0.818677 


0.818525 


0.821463 


Photon-Limited S/N (%) 


83.1 


83.1 


82.6 



''RMS frame-to-frame position difference. 
''MJD = BID - 2,450,000. 

'^We exclude frames during instrument/telescope settling, for insufficient points at a given knot, and for 
bad pixels in the photometry aperture. 
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TABLE 15 

Best-Fit Joint Eclipse Light-Curve Parameters 



P tirsmc tc r 


HD149bsll 


HD149bs21 


HD149bs31 


HD149bs32 


HD149bs33 




3.6 


4.5 


5.8 


5.8 


5.8 


rtj.i<ay ruf»iLitJii y-^., L/IA^ 


14.58 


14 57 


14.50 


14 42 


14.78 


Array Position (y, pix) 


15.66 


14.99 


14.37 


14.17 


14.67 


Position Consistency'^ i^x pix) 


0.009 


0.009 


0.020 


0.020 


0.015 


Pr\citiriTi r~"riTi CI ctpnnv ^ {^^. niv'l 


0.007 


0.004 


0.018 


0.013 


0.017 


Aperture Size (pix) 


3.75 


2.75 


2.75 


2.75 


2.75 


Sky Annulus Inner Radius (pix) 


7.0 


7.0 


7.0 


7.0 


7.0 


Sky Annulus Outer Radius (pix) 


15.0 


15.0 


15.0 


15.0 


15.0 


System Flux^ /\ ([jjy) 


528456 di 10 


334037 + 40 


216010 ± 70 


214102 it 40 


212621 it 80 


Eclipse Depth (%) 


0.040 di 0.003 


0.034 ± 0.006 


0.044 ± 0.010 


0.044 it 0.010 


0.044 it 0.010 


Brightness Temperature (K) 


2000 di 60 


1650 it 120 


1 600 ± 200 


1600 it 200 


1600 it 200 


Eclipse Midpoint'^ (orbits) 


0.5003 ± 0.0004 


0.4994 ± 0.0007 


0.502 ± 0.004 


501+'! S!!? 

"■-^ -0.005 


500+<! !'S3 


Pnlincp AAirlnninl"'^ f^lV1TT~)i tt-^^ 

l_jLllUftC iVllUUUlllL (^iVlJl^I jTr^y 


4^^^ Ri'^ft -\- n nnin 

JJ J .O / JU ZLL \J.\J\J 1 VJ 


4SQ(i Ifyfi^ -1- 001 Q 


4'^2S Q41 -1- OM 


4(^'^l -1- 010 


4001 OSQ -1-0 012 


Frlincp Mirlnninf^ ^MTD-nr-.n^ 


4515 X764 + n 0010 


4596 2676 + 001 9 


4'?25 942 + 0013 


4633 659 + 010 


4001 090 -t 01 2 

"^yyj^ .yy\j _i_ u.vjiz. 


Pr'Imcp T~liirnl"inn it a < riv<i\ 


2q + 02 




Q 2Q + 02 


Q 29 + 02 


T 29 -t 02 


Tti CTTPC c / m cttpc c T^t mp ( t-^ 1 nTC J 




n 9 -14. _|_ n n 1 7 


n 9 -14. _|_ n n 1 9 


n 7:14. -1-0 017 


n 7^4 -1-0 017 


rValillJ I^UUdllUll (^iV\^i J J 




l2l 


l2l 


I21 


I21 


IValllU, /Q 





29 ± 9 


30 it 6 


42 it 5 


26 it 5 


XvttilllJ, / Y 





7 ±4 


8 it 3 


14 it 2 


6.4 it 2.0 


JXclllllJ, 1 7 


-0 onis -1- ooos 


Q 







Q 




Yes 


Yes 


Yes 


No 


No 


Min. Number of Points Per Bin 


4 


5 


5 






Total Frames 


54080 


54080 


54080 


54080 


54080 


Rejected Frames (%) 


0.308802 


0.160873 


0.268121 


0.277367 


0.312500 


Frames Used^ 


50769 


48769 


50919 


50930 


53911 


Free Parameters 


6 


7 


7 


4 


4 


AIC Value 


50775 


48776 


155775 


155775 


155775 


BIC Value 


50828 


48837 


155924 


155924 


155924 


SDNR 


0.00280365 


0.00388070 


0.0120115 


0.0119551 


0.0119188 


Uncertainty Scaling Factor 


0.401020 


0.489317 


0.914694 


1.02508 


1.02476 


Photon-Limited S/N (%) 


93.9 


89.9 


73.6 


73.7 


74.3 



''RMS frame-to-frame position difference. 

''We multiply the HD149bs32 and HD149bs33 measured system fluxes by 0.968 to coiTect for an IRAC flux conversion issue in the S18.18 

pipeUne. 

■^Based on the period and ephemeris time given bv lKnutson eFal] {20091) . 
■"MJO = BJD - 2,450,000. 

'We exclude frames during instrument/telescope settling, for insufficient points at a given knot, and for bad pixels in the photometry aperture. 
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TABLE 16 

Best-Fit Joint Eclipse Light-Curve Parameters 



Parameter 


11 1 J 1 DS'+ 1 


rlU 1 '+VDS'f Z 


rlU 1 M-VDS'f J 


rlL^ 1 ^7 DS'+'4 


nu i^yDS J 1 


IJT~\1 /1QUcS7 
rlLJ i^yDSJZ 


Wavelength (|J-m) 


o.U 


S A 
O.U 


O.U 


H A 
O.U 


1 A 
10 


1 A 
10 


Array Position {x, pix) 


1 ^ f\fi 
i j.Uo 


14. 4U 


1 ^ A/1 
1 J.U4 


1 yl A/I 
14.04 


TA (\1 


Tl 

Zj.Vj 


Array Position (y, pix) 


14. 4j 


1 ^ AO 


14. 1 J 


1 ^ AQ 
1 J.Uo 


TT 1 A 

ZZ. 10 


TA A^ 

ZU.Oj 


Position Consistency^ i^x^ pix) 


(\ noi 


A Al 


A Al 1 
U.Ul 1 


A Al T 
U.UIZ 


A A 1 A 
U.UlO 


A Al fi 
U.Ui o 


Position Consistency^ (S\, pix) 


A (Yin 


A Al 1 
U.Ul 1 


A Al T 
U.UIZ 


A Al 1 
U.Ul 1 


A AT 1 
U.UZl 


A Al C 

U.Ulo 


Aperture Size (pix) 


A n 

4.U 


J.J 


1 7^ 
J. / J 


J.J 


7 A 
Z.U 


7 A 
Z.U 


Sky Annulus Inner Radius (pix) 


7.0 


1 A 

/.u 


1 A 

/.U 


/.U 


A A 
O.U 


A A 
O.U 


Sky Annulus Outer Radius (pix) 


1 J.U 


1 J.U 


1 J.U 


1 J.U 


i 1 .u 


iZ.U 


System Flux^ Fs (M-^y) 


11/ lyxj ^ lUU 




1 1 78 1 S -1- A 
li / olo ^ O 


1 1 eiAs -1- 70 

1 loJDJ It /U 


1 OA 1 O -L A 
loO 1 o ^ O 


1 QQAS -1- 1 1 
1 ooUJ It 1 1 


riciipse uepm y /c) 


U.UjZ it U.UUO 


A A^T -1- A AAA 
U.UjZ It U.UUO 


A A^T -U A AAA 
U.UjZ It U.UUO 


A A^T -1- A AAA 
U.UjZ It U.UUO 


A AG^ -U A AIT 
U.Uoj It U.UjZ 


A AS^ -U A ATT 
U.UoJ It U.UjZ 


Brightness Temperature (K) 


1 A^n _|_ 1 1 A 


1 A^A -I- 1 1 A 
IOjU it 1 lU 


1 A^A -1- 1 1 A 
IOjU it 1 lU 


1 A^A -I- 1 1 A 
IOjU It 1 lU 


1 GAA -1- AAA 
1 oUU it OUU 


1 SAA -1- AAA 
1 oUU It OUU 


Eclipse IVIidpoint"^ (phase) 


n sno') -1- n nnn7 

U. JUUZ ^ u.uuu / 


U.JUIU ^ \j.\)\j\jy 


n 4QA1 -1- n nni 9 

u.tyoi ^ U.UUIZ 


n AQss -1- n nnnA 

u.^+yoo ^ U.UUUO 


U.jUl J ^ U.UU1'+ 


U.jUl J ^ U.UUIM- 


liciipse iviiapoini liViJUyj^^ 


"5 AAA QAT -U A Am 

jouo.yoz It u.uuz 


A ^A7 ^ 1 "3 -1- A AAT 
4jO / . J 1 J It U.UUj 


A ^QQ 1 11 -I- A AA1 

4jyy.ijj it u.uuj 


/1Q1T AIT-LA AAT 
4yiZ.01 J It U.UUZ 


/1H7H1-I-A AA/1 

4J 1 / . J 1 1 It U.UU4 


A^AI 1 0/1 -U A AA/1 
4j4j. iy4 It U.UU4 


ticlipse Miapoint (iVlJUj^ji^) 


"JAAA OAI J- A AAO 


A ^ AT ^ 1 1 _L A AAI 

4jo/.j1j it U.UUj 


A ^QO 1 1A J- A AA1 

4jyy.ij4it u.uuj 


/1Q1T A1/1 J- A AAT 

4yiz.oi4 it U.UUZ 


/IQ 1 "7 "2 1 1 -LA AA/1 
4j 1 / .J 1 1 it U.UU4 


/IQ/1 1 1 0/1 J- A AA/1 

4j4j. iy4 it U.UU4 


Eclipse Duration (?4-] , hrs) 


1 TO J- A AT 
j.Zy it U.UZ 


1 TO _L A AT 

j.zy it U.UZ 


1 TO J- A AT 

J.zy it U.UZ 


"J TO J- A AT 

J.zy it U.UZ 


1 TO J- A AT 
J.zy it U.UZ 


1 TO J- A AT 

J.zy it U.UZ 


Ingress/Egress Time {t2-\, hrs) 


All/I 1 A CM '7 
U.ZJ4 it U.UiZ 


(\ 11 A -LA ^1 1 T 
U.Zj4 it U.UIZ 


A Tl/I 1 A A 1 T 

U.Zj4 it U.UIZ 


A TQ/I 1 A Al T 

U.Zj4 it U.UIZ 


A T1/1 J- A Al T 
U.Zj4 it U.UIZ 


A Tl/I 1 A A 1 T 

U.Zj4 it U.UIZ 


Ramp Equation {R(ty) 


f9l 


rm 


IIUI 


lzj_ 


rm 

UJJ 


f9l 


Ramp, I'o 


15.4 di 1.2 





U 


1 /I -L 
14 ± 


A 
U 


62 di 10 


Ramp, r[ 


J. i it U.J 


U 


A 

u 


A 1 A 

U it 4 


A 

u 


T/1 J- A 
Z4 it 4 


Ramp, ^2 


U 


A Ad J- A AAT 
U.USj it U.UUZ 


A AA 1 /I _L A AA 1 A 

-U.UU14 it U.UUlO 


A 
U 


A /lll-LAAIT 
U.41 J it U.UIZ 


A 

U 


Ramp, 


U 


A A7 -I- A 1 1 
-U.O / It U. 1 1 


n 
u 


A 
U 


< 1 -L A T 

-J. 1 it u.z 


U 


DLioo iviap \ivi[^x^y)) 


iNO 


les 


Woe 

ICS 


ICS 


iNO 


i>JO 


Min. Number of Points Per Bin 




Q 

o 


A 

4 


1 A 
lU 






Total Frames 


44352 


54080 


60500 


54080 


1050 


1050 


Rejected Frames (%) 


0.47574 


0.432692 


0.634711 


0.488166 


0.285714 


0.571429 


Frames Used^ 


44041 


48801 


60100 


46274 


1047 


1044 


Free Parameters 


15 


4 


3 


4 


10 


12 


AIC Value 


194243 


194243 


194243 


194245 


2113 


2113 


BIG Value 


194518 


194518 


194518 


194540 


2237 


2237 


SDNR 


0.00845740 


0.00833381 


0.00837245 


0.00847674 


0.00311603 


0.00337752 


Uncertainty Scaling Factor 


0.816657 


0.813220 


0.819215 


0.981709 


0.374578 


0.402919 


Photon-Limited S/N (%) 


82.3 


83.1 


82.4 


81.6 


27.2 


24.9 



"RMS frame-to-frame position difference. 

''We multiply the HD149bs44 measured system flux by 0.973 to correct for an IRAC flux conversion issue in the S18.18 pipeline. 
■^Based on the period and ephemeris time given by Knuts on et arU200^ . 
■^MJO = BJD - 2,450,000. 

■^We exclude frames during instrument/telescope settling, for insufficient points at a given knot, and for bad pixels in the photometry aperture. 



